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ABSTRACT 


Two-dimensional  numerical  simulation  is  a  necessary  tool  for  modern 
semiconductor  device  design.  Analytical  models  and  judicious  application 
of  one-dimensional  simulation  cannot  accurately  represent  the  highly  two- 
dimensional  impurity  profiles  and  structures  of  VLSI  devices.  Moreover,  the 
allowable  device  structures  and  bias  conditions  of  existing  two-dimensional 
simulation  programs  are  too  restrictive  to  provide  the  necessary  design  infor¬ 
mation. 

A  two-dimensional  numerical  simulation  program.  P1SCLS.  has  been 
written  in  order  to  study  various  aspects  of  device  simulation.  The  program 
uses  vectorized  LU  decomposition  to  alternately  solve  Poisson's  equation  and 
the  electron  current  continuity  equation  (Gummels  method).  The  program  is 
extremely  flexible  and  useful  in  evaluating  two-dimensional  simulation  con¬ 
cerns  such  as  grid  allocation,  boundary  conditions,  convergence  characteris¬ 
tics  and  physical  models. 

,/• 

The  discretization  grid  is  analyzed  in  comparisons  of  rectangular  and 
triangular  grids  and  in  the  allocation  of  grid  points  within  critical  regions 
of  the  device.  A  triangular  grid  achieved  by  distorting  a  rectangular  grid 
is  advocated  as  a  reasonable  compromise  between  the  flexibility  of  general 
triangular  grids  and  the  regularity  and  matrix  solution  method  computability 
of  rectangular  grids.  A  finite  difference  discretization  of  Poisson’s  equation 
and  the  current  continuity  equation  on  a  triangular  grid  is  presented. 

A  variety  of  methods  for  reducing  the  simulation  time  arc  explored.  The 
nested  dissection  grid  renumbering  scheme  is  shown  to  provide  a  significant 
storage  and  operation  count  reduction  for  larger  grids  with  a  slight  penalty 
in  vector  operation  elliciency.  Techniques  for  accelerating  the  converge  nce  of 


tv 


the  alternating  method  are  presented  which  together  reduce  solution  times 
by  a  factor  of  four  for  devices  biased  above  threshold.  These  methods 
involve  computation  of  an  improved  initial  guess,  elimination  of  excessive 
solutions  of  Poisson’s  equation,  overrelaxation  of  the  potential  updates  and 
reduction  of  the  Poisson  linearization  term.  Even  with  these  improvements, 
however,  simulations  above  threshold  still  require  about  four  times  as  long 
as  subthreshold  simulations.  This  slow  convergence  appears  to  be  correlated 
with  slow  oscillations  of  the  first  harmonic  in  spatial  frequency  of  the  surface 
potential  in  the  inverted  channel  between  source  and  drain. 

Two  application  examples  demonstrate  the  utility  of  the  PISCES  pro¬ 
gram  and  two-dimensional  numerical  simulation  in  general.  Simulation  of 
an  implanted  channel  MOSFET  reveals  a  50  fold  increase  in  punchthrough 
current  with  a  12%  increase  in  source  drain  junction  depth.  Field  depen¬ 
dent  mobility  is  investigated  with  the  implementation  of  a  distanee-from-the- 
surface  mobility  model. 
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Chapter  1 


INTRODUCTION 

As  device  geometries  shrink  in  the  pursuit  of  Very  Large  Scale  Integration 
(VLSI),  two-dimensional  numerical  simulation  of  devices  gains  in  impor¬ 
tance.  One-dimensional  approximations  which  are  valid  for  low  fields  and 
large  lateral  dimensions  with  respect  to  vertical  dimensions  no  longer  apply. 
Intensions  *  >  the  one-dimensional  theory  can  he  useful  hut  cannot  accurately 
account  for  the  highly  two-dimensional  structure  of  modern  devices. 

1.1  Perspective 


In  1070  gate  lengths  of  7  microns  were  typical  for  production  LSI  cir¬ 
cuits  [l.l].  By  I9S0  the  gate  length  had  been  reduced  to  ! 2  microns  and  cur¬ 
rent  contracts  for  the  US  Government's  Very  High  Speed  Integrated  Circuits 
(VIISIO)  program  [1.2]  call  for  production  of  0.5  micron  devices  by  19v'>. 
The  constant  reduction  of  device  dimensions  seen  over  the  last  20  years  is 
expected  to  continue  to  at  least  1900.  As  the  device  lateral  dimensions  have 
come  down  it  has  become  necessary  to  alter  other  structural  and  operational 
parameters  in  order  to  maintain  desirable  operating  characteristics.  Device 
scaling  theory  [l.d]  describes  how  In  optimally  adjust  device  vertical  dimen¬ 
sions,  voltage  levels,  and  doping  concent-rat  ions  as  a  I  unct  ion  of  lat  ri  al  dimen¬ 
sions  to  minimize  the  short  channel  effects  caused  by  strong  two-dimensional 
fields.  Practical  considerations,  however,  including  fabrication  constraints 
and  logic  level  noise1  immunity  have  caused  designers  to  sub-optimallv  scale 
these  parameters  thus  retaining  some  short-channel  behavior. 
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Figure  1.1  shows  a  comparison  of  the  equipotential  contours  of  two 
metal-oxide-semieonductor  field  effect  transistors  (MOSFET’s)  which  are 
identical  except  for  their  gate  length.  For  Figure  1.1a,  the  channel  length 
(metallurgical  junction  spacing)  is  5  microns  while  for  Figure  1.1b  it  is  2 
microns.  Note  that  in  both  cases,  as  the  equipotential  lines  near  the  surface 
curve  to  follow  the  junction  boundaries,  they  pull  away  from  the  surface  in¬ 
dicating  higher  surface  potentials  near  the  junctions  and  an  increase  in  con¬ 
trol  over  these  potentials  by  the  source  and  drain  and  a  reduction  in  control 
by  the  gate.  For  the  long-channel  device  of  (a),  these  edge  effects  are  a  small 
percentage  of  the  total  channel  length  and  thus  have  limited  influence  on 
device  characteristics.  In  the  short-channel  device,  however,  the  edge  effects 
extend  throughout  a  large  percentage  of  the  channel  length  and  have  a  strong 
influence  on  the  device  characteristics.  One  result  of  this  is  a  lower  threshold 
voltage  than  that  predicted  by  theory  since  the  surface  potential  is  higher 
for  a  given  gate  bias.  Another  result  is  an  increased  sensitivity  of  the  output 
current  to  the  drain  bias  (drain  conductance)  due  to  the  control  exerted  by 
the  drain  on  the  surface  potential  in  the  channel. 

Such  effects  are  clearly  the  result  of  the  two-dimensional  structure,  and 
attempts  have  been  made  to  model  these  ('fleets  analytically.  Some  success 
has  been  obtained  by  Yau  [l.-lj  and  extended  by  others  [1.5]  to  model  short- 
channel  effects  by  using  a  charge  sharing  theory.  The  basic  features  of  this 
model  are  shown  in  Figure  1.2.  The  charge  controlled  by  tin'  gate  is  assumed 
to  be  contained  within  the  trapezoidal  region  immediately  beneath  the  gate 
with  the  remaining  charge  controlled  by  the  source  and  drain  respectively. 
Further,  the  junctions  and  depiction  edge  boundaries  are  assumed  to  he 
cylindrical. 
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ng.  1.1.  A  comparison  of  short-  and  long-channel  MOSFKT  oquipot  out  ial 
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This  model  roughly  approximates  (ho  threshold  shifts  seen  a-,  tin1  chan¬ 
nel  length  is  reduced  in  devices  fabricated  using  conventional  long-channel 
techniques;  however,  as  device  structures  are  optimized  for  short  channel 
performance,  the  assumptions  made  (o.g.  cylindrical  junctions)  no  longer 
apply  and  the  model  is  invalidated.  Thus  with  analytical  modeling,  if  all  of 


> 

i 


3 


GATE 

VZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZk 


r»— oxide 


SOURCE 


\  GATE  CHARGE  // 


DRAIN 


>- 


■-< 


SOURCE  CHARGE  / 


/ 


\ 


\ 


DRAIN  CHARGE 


DEPLETION  EDGE 


SUBSTRATE 


Fig.  1.2.  Charge  sharing  modal  of  Yau. 

the  assumptions  made  are  valid,  then  the  proper  results  will  he  obtained; 
however,  if  the  assumptions  break  down  then  the  results  will  be  invalid.  The 
principal  advantage  of  numerical  simulation  of  semiconductor  devices  is  that 
few  a;  umptions  need  to  be  made.  If  there  are  design  errors  which  result 
in  excessively  large  electric  fields  or  alternate  conduction  paths,  for  example, 
the  simulation  will  accurately  reveal  the  existence  of  thee  ('fleets  even  though 
they  had  not  been  expected.  The  principal  disadvant ages  of  numerical  simula¬ 
tion  are  that  there  are  no  simple  equations  to  describe  device  behavior,  it  is 
not  always  obvious  how  to  interpret  the  results,  and  a  significant  amount  of 
computation  is  required. 

or  course,  attempts  can  and  will  be  made  to  model  new  structures. 
Analytical  models  provide  insight  into  the  interplay  of  device  parameters 
which  can  be  useful  to  device  designers  and  users;  however,  the  modeling  job 
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becomes  increasingly  dillicult  and  the  assumptions  become  more  restrictive  as 
device  structures  shrink  in  size  and  become  more'  complex.  In  fact,  numerical 
simulation  is  inevitably  used  in  arriving  at  or  verifying  analytical  models. 

As  an  example  of  the  complexity  of  modern  device  structures,  Figure  1.3 
shows  two  exotic  field  ellect  devices:  a  static  induction  transistor  (l.Oj  and 
a  taper-isolated  dynamic-gain  RAM  cel!  [1.7].  The  operating  characteristics 
of  these  devices  result  primarily  from  the  two-dimensional  nature  of  their 
impurity  profiles  and  electric  fields.  The  creation  of  analytical  models  for 
these  devices,  applicable  throughout  their  entire  operating  regime,  would 
clearly  be  a  difficult  task. 

For  sonic  applications,  (e.g.  circuit  simulation  programs)  empirical  rela¬ 
tions  suffice  for  explaining-  device  behavior.  This  form  of  modeling  suffers 
from  several  drawbacks,  f  irst,  the  model  parameters  often  have  no  physical 
basis  and  may  only  be  extracted  from  measured  data.  Second,  since  each 
region  of  device  operation  requires  a  different  empirical  relation,  it  is  dillicult 
if  not  impossible  to  match  the  device  character;:  ! ics  in  the  t ran -it ion  from  one 
region  to  another.  Going  a  step  further,  one  recent  circuit  simulation  pro¬ 
gram  uses  tallies  instead  of  closed  form  expression.-*  for  device  charact crist ics 
[I.S],  Much  physical  basis  is  lost  in  this  method  except  in  the  determination 
of  how  to  parameterize  the  tables.  In  either  ease,  closed-form  expressions 
or  tables,  numerical  simulation  can  be  used  to  generate  the  required  device 
characteristics.  Numerical  simulation  also  provide-,  a  significant  opportunity 
to  study  the  edi  cts  of  dc\ ice  technology  variables  on  circuit  performance. 
Furthermore,  the  use  of  process  simulation  to  generate  the  device  profiles 
provides  the  opportunity  to  directly  study  the  link  between  fabrication  steps 
and  circuit  performance. 
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fed  back  into  the  process  specifications  and  the  sequence  is  repeated  until 
satisfactory  agreement  is  obtained.  The  device  lot  wafers  are  then  processed 
and  electrical  measurements  are  made  from  which  two-dimensional  profile 
information  may  be  inferred.  Electrical  measurements  are  also  made  to 
determine  device  and  circuit  characteristics.  Feedback  is  provided  at  each 
level  to  allow  optimization  of  the  process. 

The  principal  savings  in  cost  and  time  for  simulation  versus  actual 
fabrication  and  testing  of  devices  comes  in  process  simulation.  Typically, 
one  simulation  of  a  process  using  Sl'FllHM  would  take  only  a  few  minutes, 
costing  tens-of-dollars  on  a  main  vamo  compute'!'.  Actual  fabrication  would 
typically  take  several  weeks  and  co.-t  thousands  of  dollars.  Obviously,  the 
savings  re  idling  from  the  use  of  -iinnlation  are  substantial.  The  situation 
is  reversed  somewhat  for  device  and  circuit  simulation  verm-  device  and  cir¬ 
cuit  measurements.  Device  and  circuit  electrical  measurements  are  relat i\ ely 
quick  and  inexpensive  ta-  l;s  while  the  cost  of  each  simulation  is  roughly  com¬ 
parable  to  that  for  process  simulation.  Even  so,  device  simulation  offers  a 
tremendous  advantage  over  device  measurement  since  one  him  <1  noi  I’nbri.  ate 
the  device  first.  Moreover,  simulation  provides  a  detailed  two-duwn-ional 
view  of  the  physics  which  device  measurements  cannot  provide. 

Another  advarit  age  oT  si  mu  lat  ion  lies  in  1  lie  fact  that  it  is  oft  cn  di  I  limit  1  o 
accurately  measure  two-dimensional  device  structures  [i.l  Ij.  When  attempt¬ 
ing  to  analyze  or  improve  device  performance,  one  would  like  to  measure 
device  profiles.  In  case.;  where  accurate  measurements  are  not  possible, 
simulation  allows  the  engineer  the  opportunity  to  manipulate  the  simulated 
profiles  and  observe  the  effects  on  device  charnel  eristics,  thus  inferring  the 
actual  profile  shapes  [l.l‘2|. 

Finally,  simulation  provides  a  much  greater  insight  into  device  behavior 
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Fig-  1.1.  The  parallel  paths  of  simulation  versus  actual  wafer  proceed ng. 

Feedback  is  provided  at  each  level  by  the  comparison  of  measured 
or  simulated  result:-  on  profiles  or  operating  characteristic.--  against 
the  desired  values. 


than  can  he  obtained  in  any  other  way.  In  the  case  of  unwanted  punch!  hrough 
currents  in  MOSFFTs  [1.13],  for  example,  one  can  both  determine  terminal 
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conditions  of  punchthrough  and  spatially  identify  the  point  of  punchthrough. 
Further,  a  proposed  solution  to  the  problem  may  he  simulated  arid  its  success 
or  failure  determined. 

1.2  History 

The  history  of  numerical  simulation  of  semiconductor  devices  begins 
with  the  work  of  Gumme!  [l.l-l]  in  19t>l  on  the  ono-dimen-ional  steady-state 
analysis  of  bipolar  transistors.  Ills  method  provided  a  two-carrier  solution 
meaning  that  the  continuity  and  transport  equations  were  satisfied  for  both 
holes  and  electrons.  The  principal  contribution  of  his  early  work  was  the  use 
of  a  sequent  ini  iteration  sc  In' me  for  obt  aining  a  consist  out  so!  ut  ion  to  the  three 
sets  of  equations:  Poisson’s  equation,  and  the  h'le  and  electron  continuity 
equations,  hi  his  method,  one  first  solves  Poisson's  equation  followed  by 
the  electron  and  hole  continuity  equations  in  succession.  The  cycle  is  then 
repeated.  The  solution  of  each  equation  individually  in  this  manner  requires 
much  less  work  t  Inn  eo!\  ing'  all  t  liree  simultaneously.  Generally,  however,  t  he 
convergence  is  not  as  fast  as  the  quadratic  convergence  which,  can  be  obtained 
for  a  simultaneous  solution  [1.15].  Nevertheless,  the  amount  of  work  saved 
in  each  iteration  generally  compensates  for  the  slower  convergence.  These 
points  will  be  examined  in  greater  detail  in  Chapter  1. 

l)e  Mari  analyzed  the  p-n  junction  in  one  dimension  in  19(1S  [1.1(1  and 
enhanced  the  numerical  analysis  capabilities  to  include  transient  conditions 
[1.17].  In  1999  SchaiTrllor  and  Gummel  published  their  work  on  the  transient 
analysis  of  a  Read  diode  oscillator  [I.ISj.  This  paper  provided  the  second 
major  advance  in  numerical  algorithms  with  the  introduction  oT  a  carrier 
transport  equation  descrot i/.at ion  scheme  which  allowed  larger  grid  -pacing 
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and  llms  fewer  variables.  This  method  will  be  described  further  in  Chapter 
3.  The  advent  of  two-dimensional  simulation  in  1969  lowered  interest  in  one¬ 
dimensional  simulation,  thus  slowing  its  development  and  focusing  its  ap¬ 
plication  to  bipolar  devices.  The  most  significant  works  to  follow  emphasized 
specific  device  analyses  obtained  from  one-dimensional  simulations  and  the 
inclusion  of  higher  order  physical  phenomena  such  as  band  gap  narrowing  and 
mobility  variations  [1.19,  1.20].  Selected  woi  ks  in  one-dimensional  analysis 
include  Gokhale  in  1970  [1.21];  llachtel,  et  al.  in  1972  [l.22;;  and  D'Avanvo 
in  1979  [1.23]. 

Two-dimensional  simulation  appeared  in  the  literature  in  19fi9  with  the 
publication  of  works  by  Kennedy  and  O'Brien  1.21  20  on  the  simulation  of 
JI'fiTs  and  by  Slotboom  [1.27,  1.28]  on  bipolar  ttan-i  tors.  The  advent  of 
short  channel  IT'.Ts  in  this  period  was  the  driving  factor  for  two-dimen  ion.d 
simulation  and  almost  all  subsequent  work  was  aimed  at  ltd  i.Ts.  Rei-er 
introduced  two-dimensional  transient  anly-ds  in  i<)7()  [1.29  3  l]  and  Mock 
presented  his  stream  function  formulation  of  the  earlier  transport  equations 
in  1973  [1.31]. 

The  next  dgnificant  development  came  in  1973  with  the  first  publication 
of  finite  element  analysis  by  llachtel  jl.3f>,  1.3(5]  followed  in  197  1  by  Barnes 
and  Lomax  [1.37,  1.38]  and  Butorin  and  Cottrell  1.39  11].  \  1 1  previous 

work  had  been  based  on  finite-di ITerenee  discretization  schemes  ami  their  as¬ 
sociated  rectangular  grids.  ’The  use  or  finite  elements  pmCuled  improvement- 
over  finitc-dilTcronce  discretization  with  the  ability  to  model  non- n  cl  angular 
structures  and  a  more  efficient  use  of  grid.  Although  rectangular  grids  do 
not  prohibit  the  simulation  of  noil-rectangular  structures  [1.12],  the  numeri¬ 
cal  techniques  for  accomodating  these  structures  had  not  been  implemented. 
Hence,  previous  work  had  been  restricted  to  reel  angular  structures.  These 


topics  will  be  discussed  further  in  Chapter  2. 

Developments  since  1975  have  focused  on  decreasing  program  size  and 
solution  time  while  increasing  the  accuracy  of  the  solution  in  terms  of  both 
the  numerical  algorithms  used  and  the  models  of  the  device  physics.  Another 
aim  has  been  the  development  of  “friendly”  user  interfaces  for  the  simulation 
programs  in  order  to  make  the  simulation  technology  more  available  to  device 
designers  [1.13].  This  progression  from  use  of  simulation  as  a  laboratory 
tool  in  developing  semiconductor  device  theory  to  application  in  produc¬ 
tion  facilities  for  optimizing  device  structures  is  a  field  still  in  its  infancy. 
Nonetheless,  its  application  holds  great  promise  in  providing  the  increased 
productivity  needed  for  VLSI  design. 

Several  programs  have  recently  become  widely  available.  NliMOS,  based 
on  the  early  work  of  Kennedy  [  1 .21];  CADDKT,  based  on  Mock's  work  [l.-Hj; 
and  MINIMOS,  [1.15]  all  provide  steady-state  solutions  for  essentially  rec¬ 
tangular  l'KT  structures.  TWIST  [l.-lfij  and  ClMMlXl  [1.-12]  solve  only  for 
the  device  potentials  but  are  quite  useful  for  simulation  of  sub-threshold 
and  punchthrough  characteristics  as  well  as  device  breakdown  characteris¬ 
tics.  TWIST  is  limited  to  rectangular  geometries  while  CLMIXI  allows  lion- 
rectangular  structures. 

Currently,  work  is  in  progress  on  two  fronts  which  will  provide  valuable 
aid  to  the  device  designer  the  development  of  three-dimensional  simula¬ 
tion  and  the  the  introduction  of  simplified  fast  two-dimensional  simulation. 
Several  researchers  have  published  preliminary  work  on  three-dimensional 
simulation  [1.47-49].  Most  notable  of  these  is  that  of  Buturla  and  Cottrell 
with  their  extension  of  the  two  dimensional  simulation  program  1  ll.LDAY 
[1.50]  to  three  dimensions.  These  simulations  have  shown  that  there  are 
characteristics  of  short  and  narrow  semiconductor  devices  which  can  he  ac- 


curat ely  simulated  only  in  three  dimensions.  These  programs,  however,  must 
run  on  largo  mainframe  computers  and  consume  considerable  computer  time; 
thus,  they  are  currently  of  limited  practical  value  tc  device  designers. 

At  the  other  end  of  the  spectrum  is  the  simplified  two-dimensional 
program  of  Oh,  SDV1CE  [1.51].  This  program  provides  a  very  fast  solu¬ 
tion  to  the  two-dimensional  FET  transient  problem  and  uses  very  little 
computer  memory.  It  is  limited  in  the  device  geometries  and  operating 
regions  which  it  handles  and  the  accuracy  of  its  solutions  require  further 
verification.  However,  this  program  solves  a  particular  type  of  problem  ex¬ 
tremely  efficiently. 


1.3  PISCES 


A  two-dimensional  numerical  simulation  program  has  boon  written  for 
the  purpose  of  investigating  grid  and  boundary  condition  sens'll ivites,  con¬ 
vergence  limitations,  device  physics  models,  and  to  compare  other  numeri¬ 
cal  simulation  programs.  PISCES  (Poisson  and  Single-carrier  Continuity 
Equation  Solver)  solves  the  Poisson  equation  and  the  steady-state  electron 
continuity  equation  using  the  alternating  method  (Ciunimel’s  algorithm)  on  an 
IIP-  1000P  minicomputer.  The  program  uses  a  finite  difference  discretization 
on  an  irregular  triangular  grid  and  thus  easily  handles  non  planar  surfaces 
and  interfaces. 

The  program  was  written  for  use  on  field  effect  transistors  where  a 
single-carrier  solution  is  sufficient.  Field  effect,  transistors  are  majority  carrier 
devices  and  very  little  error  is  introduced  by  ignoring  the  minority  carriers 
except  in  extreme  biasing  conditions  such  as  avalanche  breakdown.  Hy  solv¬ 
ing  only  the  electron  continuity  equation,  time  is  saved  since  the  hole  con- 
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tinuity  equation  does  not  have  to  be  solved,  and  slightly  faster  convergence 
is  obtained. 


The  equations  solved  are: 


Poisson 
Continuity 
Transport 
Ilolt  7,  man  n 


V  •  (eVVO  =  n  —  p  —  AT 
V  •  Jn  =  q(Rn  -  Gn) 

Jn  =  -qnnTiVip  +  qDnVn 
n  =  nieq^-^'>/kT 


p  ~  n^-M/kT 


where  (  is  tlu’  permittivity,  ip  is  the  electric  potential,  n  and  p  are  the  free 
electron  and  hole  concentrations,  iV  is  the  net  ionized  impurity  concentre!  ion, 
J n  is  the  electron  current  density,  Rn  and  Gn  are  the  electron  recombination 
and  generation  rates,  /<„  is  the  electron  mobility.  Dn  is  the  electron  diffusion 
constant,  ?i,  is  the  intrinsic  carrier  concentration,  and  6V  are  the  eledron 
and  hole  quasi-Fermi  levels,  and  kT  /  q  is  the  thermal  voltage. 

The  extension  of  finite  difference  methods  to  triangular  grids  is  a  novel 
approach  for  semiconductor  device  simulation  although,  it  has  been  applied 
in  other  fields  [1.52].  Irregular  triangular  grids  possess  desirable  features 
for  semiconductor  simulation  including  the  ability  to  conform  to  irregular 
shapes  and  the  capacity  for  local  grid  refinement  without  inducing  excessive 
grid  elsewhere.  These  same  advantages  are  the  driving  factors  which  have 
resulted  iu  the  development  of  finite'  element  simulators. 

The  IMSCIvS  program  also  contains  user  oriented  features  which  increase 
its  utility.  I  sing  the  input  parser  and  graphics  interface  routines  developed 
for  the  (IKMINI  program,  P1SCKS  provides  a  flexible,  friendly  user  interface 
for  both  input  and  output.  Two  examples  of  contour  plot  output  from 
PfSCFS  are  shown  in  f  igure  1.5.  ICquipotcril ial  contours  arc  shown  for  a 


typical  short-channel  MOSFPT  and  for  a  MESFICT.  A  complete  PI  SC  ICS 
simulation  example  is  provided  in  Appendix  A. 

1.4  Overview 

The  objective  ot'this  work  is  to  aid  the  development  of  device  simulation 
by  examining  various  aspects  of  the  problem.  A  versatile  two-dimensional 
simulation  program  has  been  developed  in  the  course  of  this  study.  The  ver¬ 
sa!  il’ty  is  achieved  through  choice  of  grid  and  discretization  schemes  which 
allow  simulation  of  modern,  highly  two-dimensional,  non-planar  semiconduc¬ 
tor  devices  on  a  minicomputer.  The  use  of  a  minimum  number  of  grid  points 
and  application  of  several  novel  convergence  acceleration  techniques  reduces 
solution  times  to  practical  limits. 

The  effects  of  grid  on  two  dimensional  numerical  simulation  of  semicon¬ 
ductor  devices  are  discussed  in  Chapter  2.  The  various  types  of  grid  are 
presented  and  their  applicability  to  device  structures  and  solution  methods 
are  considered.  The  tradeoff  of  number  of  grid  points  versus  complexity 
of  solution  method  is  addressed.  Requirements  for  high  grid  density  in  lo¬ 
calized  regions  are  also  considered.  The  chapter  concludes  with  a  discus¬ 
sion  of  boundary  condition  sensitivites  and  the  requirements  for  accurate 
discretization  of  impurity  profiles  including  lateral  extensions  of  source  and 
drain  regions  in  the  simulation  window. 

Finite  difference  discretization  on  an  irregular  triangular  grid  is  the  sub¬ 
ject  of  Chapter  .'5.  The  discretization  and  linearization  of  Poisson's  equation 
with  carrier  statistics  constraints  is  derived  including  the  special  case  of  dis¬ 
cretization  when  tin-  grid  contains  obtuse  triangles.  Discretization  of  the 
electron  continuity  equation  is  also  described.  The  uoti-exPl  em-e  of  an  ex- 
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act  two-dimensional  equivalent  to  the  one-dimensional  Gummel-Scharfet ter 
scheme  is  proved  and  a  quasi-two-dimensional  discretization  is  described. 

Chapter  4  addresses  methods  for  solving  the  discretized  equations. 
Techniques  such  as  the  Fast  Fourier  Transform,  the  conjugate  gradient 
method,  and  relaxation  methods  are  discussed  in  connection  with  solving  the 
matrix  equations  resulting  from  the  discretization  of  Poisson’s  equation  or 
the  continuity  equation.  Methods  for  reducing  the  amount  of  required  com¬ 
putation  by  renumbering  the  grid  are  examined.  The  various  ways  of  solving 
the  set  of  coupled  equations  are  described  with  emphasis  on  the  alternating 
method  and  its  convergence  properties.  The  convergence  rate  is  shown  to  vary 
with  device  operating  conditions  and  with  carrier  mobilty.  Several  methods 
of  accelerating  the  convergence  of  the  alternating  met  hod  are  presented. 

Chapter  5  presents  two  program  application  examples.  The  first,  is  a 
typical  application  in  device  design  in  which  punchthrough  current  is  shown 
to  vary  greatly  with  a  change  in  source/drain  junction  depth.  The  second  ap¬ 
plication  concerns  mobility  phenomena  observed  in  strong  inversion.  Both  of 
these  applications  demonstrate  two  important  benefits  of  numerical  analysis 
programs  for  semiconductor  device  design.  First,  these  programs  can  be  an 
aid  in  developing  and/or  proving  theories  about  device  behavior.  Second,  no 
a  priori  knowledge  of  device  operating  conditions  are  required.  This  is  in  con¬ 
trast  to  the  analytical  modeling  case  where  the  proper  analytical  model  must 
be  chosen  depending  on  the  device  region  of  operation  (e.g.  subthreshold, 
linear,  breakdown). 

The  conclusions  of  this  research  and  recommendations  for  further  work 


are  contained  in  Chapter  6. 


Chapter  2 


GRID 

The  execution  time  and  storage  requirements  of  two-dimensional  semi¬ 
conductor  device  simulation  programs  are  directly  dependent  on  the  number 
of  grid  points  (nodes)  in  the  discretized  analysis  space.  The  number  of  equa¬ 
tions  to  be  solved  is  generally  linearly  related  to  the  number  of  nodes  and  the 
number  of  arithmetic  operations  required  for  the  solution  is  proportional  to 
Nn  where  N  is  the  number  of  nodes  and  a  is  somewhere  between  1.5  and  2. 
Reducing  the  number  of  nodes  in  a  simulation  is,  therefore,  a  matter  of  great 
importance. 

2.1  Grid  Types 

There  arc  two  types  of  grid  which  are  of  interest  in  two-dimensional 
device  simulation:  rectangular  and  triangular.  Within  each  typo  there  are 
variations  which  have  substantial  impact  on  the  number  of  nodes  and  on 
the  solution  methods  which  can  be  used.  Figure  2.1  shows  several  of  the 
principal  variations.  Figure  2.1a  shows  a  regular  rectangular  grid  in  which 
the  grid  spacing  is  constant,  although  not  necessarily  the  same,  in  both  the 
vertical  and  horizontal  directions.  This  grid  has  the  desirable  feature'  that 
the  discretization  coefficients  are  constant  in  both  directions  so  storage  is 
minimized.  It  is  also  the  grid  on  which  nearly  all  numerical  analysis  theory 
is  based,  thus  it  was  the  grid  used  for  most  of  the  early  work  on  device 
simulation.  This  grid  is  suitable  for  any  of  the  solution  methods  discussed  in 
the  next  chapter.  Unfortunately,  the  semiconductor  device  problem  requires 
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that  the  grid  be  very  fine  in  some  regions  of  the  device,  but  not  necessarily 
so  in  others.  Therefore,  the  constant  spacing  rectangular  grid  wastes  a  lot  of 
nodes  in  regions  of  the  device  where  fine  spacing  is  not  required. 

The  semi-constant  spacing  rectangular  grid  of  Figure  2.1b  partially 
resolves  this  problem  by  allowing  the  grid  to  have  variable  spacing  in  one 
dimension.  The  savings  in  grid  are  not  substantial,  however,  so  this  grid  is 
of  no  great  consequence  except  in  its  applicability  to  Fast  Fourier  Transform 
(FFT)  solution  techniques.  The  uniform  grid  spacing  in  the  horizontal  direc¬ 
tion  supports  a  fast  solution  to  Poisson’s  equation  through  the  use  of  the  FFT. 
This  technique,  its  advantages  and  limitations,  will  be  discussed  further  in 
Chapter  4. 

The  most  common  grid  is  the  variable  spacing  rectangular  grid  of  Figure 
2.1c.  The  grid  spacing  is  allowed  t.o  vary  in  both  directions,  yet  the  coefficient 
storage  required  for  an  m  by  n  grid  is  only  on  the  order  of  w  +  n.  The 
uniformity  of  the  overall  structure  allows  for  simple,  straight-forward,  easy- 
to-program  algorithms  for  equation  solution  regardless  of  the  solution  method 
used.  The  sole  disadvantage  of  this  grid  is  that  it  is  still  rather  inefficient  in 
grid  allocation.  Fine  grid  spacing  at  any  point  within  the  device  results  in 
grid  lines  which  extend  this  spacing  throughout  the  device  in  horizontal  or 
vertical  bands. 

The  grid  of  Figure  2. Id  is  a  special  case  of  a  class  of  grids  in  which  the 
grid  lines  are  terminated  within  the  simulation  region  of  the  device.  The 
uniform  horizontal  spacing  of  this  grid  and  the  fact  that  the  spacing  remains 
uniform  and  exactly  doubles  as  the  vertical  grid  lines  terminate  means  that 
this  grid  is  also  aplicablc  to  FFT  solution  techniques.  A  more  generalized 
form  of  this  grid  has  variable  spacing  in  both  directions  and  grid  lines  which 
may  terminate  in  either  direction.  This  is  the  type  of  grid  used  by  Adler  [2.1) 
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i/i  his  simulation  of  thyristors.  This  grid  is  obviously  flexible  in  its  ability  to 
place  regions  of  coarse  and  fine  grid  throughout  the  device.  The  principal 
shortcoming  of  this  type  of  grid  is  that  as  the  number  of  terminating  grid 
lines  increases,  most  of  the  advantage  of  rectangular  grids  (i.c.  small  storage 
anil  simple  algorithms)  is  lost  and  one  may  as  well  use  a  triangular  grid. 

The  triangular  grid  of  Figure  2.1e  is  based  on  a  rectangular  grid  which 
has  been  distorted  so  that  it  conforms  with  features  of  the  device  diagonals 
are  added  to  divide  each  rectangle  into  two  triangles.  This  grid  has  several 
desirable  features.  First,  the  underlying  rectangular  grid  is  easy  for  a  user 
to  specify.  Moreover,  it  is  not  difficult  to  specify  operations  which  distort 
the  grid  to  the  desired  shape.  Second,  this  grid  retains  its  rectangular 
connectivity  and  thus  supports  simple  solution  method'  such  as  line  iterative 
techniques  and  maintains  a  well- ih  lined  matrix  structure.  On  the  negative 
side,  each  node  has  a  unique  set  of  discretization  coellicients  so  that  storage 
for  an  in  by  n  grid  is  on  the  order  of  mn,  the  number  of  nodes.  Another 
dilliculty  with  this  grid  is  that  in  regions  of  groat  distortion,  the  triangles  may 
become  unavoidably  obtuse.  This  condition  should  he  a'.oided  if  pe.  able. 

Finally,  the  completely  genera!  triangular  grid  of  Figure  2.  If  is  the  most 
efficient  grid  allocation  in  number  of  nodes.  Spacing  may  be  made  arbitrarily 
fine  or  coarse  in  any  local  region  with  no  global  impact  except  that  the 
transition.,  between  region.1;  should  be  gradual,  (’oeflicient  storage  is  again  on 
the  order  of  the  number  of  nodes  mn.  'The  regular  structure  of  all  previous 
grids  is  lost,  however,  so  that  the  ..olution  techniques  are  nccc .-mily  less 
structured.  The  I’ISCFS  program  was  written  assuming  this  type  of  grid: 
however,  the  grid  generation  portion  of  the  program  generates  the  more 
structured  grid  of  Figure  2.1e. 

The  discussion  of  grid  types  to  this  point  has  as-aimed  simulation  of 
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Fig.  2.2.  Ov erlnying  grid  on  non-root  angular  structures. 

a  rectangular  struct  me-  with  regions  of  varying  internal  grid  density.  A 
more  realistic  case  for  semiconductor  device  simulation  is  to  allow  the  device 
structure  to  he  non-rect  angular.  This  is  mandatory  for  the  modern  two- 
dimensional  devices  described  in  Chapter  1.  Figure  2.2  shows  four  me!  lied - 
of  overlaying  a  grid  on  a  non-rect  angular  structure.  The  structure  shown  is 
the  source  region  of  an  lCl'FT  with  its  overlying  tapered  oxide. 

In  Figure  2.2a  the  taper  in  the  oxide  is  simply  ignored  and  the  device 
is  assumed  rectangular.  The  implication  is  that  the  effect  of  the  portion  of 
the  device  outsioe  of  the  simulation  region  is  insignificant.  '1'his  generally 


is  not  a  good  assumption  although  nearly  all  simulation  programs  using  the 
finite  difference  method  have  used  this  type  of  simulation  region.  Making 
this  assumption  also  requires  that  the  insulator  be  modeled  as  a  planar  slab. 
Programs  such  as  CADDET  have  taken  advantage  of  this  with  the  further 
as&uruption  that  the  electric  field  in  the  insulator  is  one-dimensional.  These 
assumptions  simplify  the  simulation  process  and  allow  a  more  rapid,  albeit 
less  accurate,  solution.  Also,  the  FFT  method  mentioned  earlier  demands 
this  type  of  grid  overlay  since  it  requrios  that  the  grid  must  be  rectangular 
with  uniform  horizontal  spacing.  In  addition,  certain  physical  parameters 
such  as  permittivity  must  he  constant  along  any  horizontal  grid  line,  limiting 
all  material  interfaces  to  he  planar  surfaces. 


Figure  2.21)  shows  an  overlay  alternative  in  which  the  rectangular  grid 
structure  is  maintained  at  the  expense  of  wasting  nod-'s  which  He  e\l<iual 
to  the  device  (i.e.  in  the  space  above  the  thin  gate  oxide).  This  appears 
to  be  a  useful  approach  when  using  solution  algorithms  which  depend  on 
a  rectangular  simulation  region.  No  application  of  this  type  of  grid  has 
been  found  in  the  literature.  When  the  solution  algorithms  do  not  require  a 
rectangular  solution  region,  t In'  nodes  externa!  to  the  device  may  be  ignored, 
and  tin'  grid  of  Figure  2.2c  results.  This  is  the  type  oT  grid  ue  d  by  the 
GEMINI  simulation  program.  Note  that  both  grids  (h)  and  (c)  require  an 
increased  grid  density  at  curved  boundaries  in  order  to  accurately  represent 
the  boundary  shape.  Finally,  the  fully  conforming  triangular  grid  of  Figure 
2. 2d  provides  the  most  accurate  simulation  with  the  minimum  number  of 
nodes. 


2.2  Criil  Density  Criteria. 


Given  a  flexible  discretization  i <i ,  one  must  then  have  a  criteria  for 
the  proper  placement  of  the  grid  points  within  the  simulation  region.  In 
semiconductor  device  simulation,  there  are  two  principal  driving  factors: 
accurate  representation  of  the  potential  and  accurate  representation  of  the 
net  charge. 

In  the  discretization  of  Poisson’s  equation,  the  assumption  is  made  that 
the  potential  varies  linearly  between  nodes  (i.c.  (In'  electric  field  is  con-taut), 
thus  the  grid  spacing  must  be  made  sufficiently  small  in  each  direction  so 
that  a  piecewise  linear  approximation  to  the  true  continuous  potential  is 
sufficiently  accurate.  This  implies  that  the  grid  must  be  the  finest  in  regions 
of  high  curvature  of  the  potential.  Prom  Poisson’s  equation,  it  i->  easily  seen 
that  regions  of  high  curvature  of  the  potential  correspond  to  regions  of  high 
net  charge  density, 
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lor  IGI'TCTs,  the  net  charge,  may  he  large  in  the  surface  inversion  layer 
where  there  are  large  numbers  of  free  carriers  or  in  depletion  regions  where 
ionized  impurities  dominate.  Vertical  grid  spacing  in  the  in\cr.-ion  layer  is 
typically  .01  //m  or  less  at  the  surface.  Neutral  regions,  no  mailer  how  highly 
doped  with  impurities,  do  not  generally  requin'  dense  grid  concent  nations 
unless  there  is  some  chance  that  charge  depletion  or  accumulation  may  occur 
there.  The  adequacy  of  a  particular  grid  spacing  for  a  given  problem  may 
be  determined  qualitatively  by  plotting  the  potential  versus  distance  in  the 
device  and  observing  whether  or  not.  the  piecewise  nature  of  the  potential 
approximation  is  evident.  The  same  check  may  he  performed  quantitatively 
by  comparing  first  and  second  order  polynomial  curve  (its  to  the  discrete 
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potential  values. 


It  is  also  desirable  that  the  total  potential  change  between  nodes  not 
be  too  large.  The  criteria  for  potential  changes  between  nodes  depends  on 
the  device  characteristics  being  simulated.  For  example,  Figure  2.3  shows 
potential  plotted  versus  lateral  distance  along  the  insulator-semiconductor 
surface  for  a  device  biased  in  punchthrough.  The  source  is  at  the  left,  the 
channel  region  in  the  middle,  and  the  drain  at  the  right.  For  this  bias 
condition,  there  is  no  inversion  layer  and  conduction  is  impeded  by  the 
potential  barrier  near  mid-channel.  The  barrier  height  is  strongly  controlled 
by  the  gate;  however,  it  is  also  under  control  of  the  drain  by  virtue  of  the 
large  drain  bias.  The  drain  current  is  proportional  to  eVa"T  where  \'jj  is 
the  barrier  height  and  Vr  is  the  thermal  voltage  (npproximtely  2(5  millivolts). 
Thus,  a  10  millivolt  error  in  the  simulated  barrier  height  can  result  in  an 
error  in  the  simulated  punchthrough  current  of  nearly  fi()f  c depending  on  grid 
spacing  and  device  profiles.  As  can  be  seen  in  llie  figure,  the'  potential  drops 
approximately  3  volts  in  .3  pm  near  the  drain;  therefore,  the  horizontal  grid 
spacing  in  this  region  must  be  fine  enough,  that  less  than  10  millivolts  error 
is  made  in  discretizing  the  3  volt  drop.  An  error  of  30',  in  subthreshold 
and  punchthrough  currents  is  typical  for  simulation  programs  due  to  this 
sensitivity.  A  similar  situation  exists  for  av.dance  breakdown  simulations 
where  large  voltage  drops  exist  across  short  dh  lances.  In  this  case,  however, 
one  is  generally  most  int  crest  ed  in  I  lie  volt  age  at  which  breakdown  occurs,  not 
in  accurate  simulation  of  the  current  near  breakdown.  The  voltage  accuracy 
required  is  also  not  critical,  being  on  the  order  of  a  few  percent .  The  voltage 
drop  per  grid  space  may,  therefore,  be  fairly  large.  Fortunately,  currents 
are  not  as  sensitive  to  voltage  errors  in  the  linear  and  saturation  regions  of 
operation  where  accurate  current  estimates  are  normally  required. 
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Hg.  2.3.  Potential  along  the  insulator-semiconductor  interface  for  a  device 
in  punr.hthrough  (VG  —  -.ll7,  Vna  —  3K,  Vs  —  OK,  V[}  =-  IK). 

If  the  region  of  high  charge  density  is  very  thin  then  the  change  in 
potential  across  that  region  will  he  small  in  spile  of  the  large  curvature 
since  potential  is  proportional  to  the  second  integral  of  charge.  In  this  ease, 
the  requirement  for  accurate  representation  of  net  charge  dictates  the  grid 
density  required.  This  is  typical  of  surface  inversion  layers  where  mobile 
charge  densities  may  vary  l>v  several  orders  of  magnitude  over  a  distance  of 
.1 //tn  or  less.  Since  KjFKT  drain  current  is  proportional  to  the  integral  of 
channel  charge  vertically  from  the  insulator  to  the  neutral  bulk,  it  follows 
that  one  needs  a  liner  grid  in  regions  of  high  net  charge  concentration  than 
in  regions  of  low  net  charge.  Thai,  is,  a  10%  error  in  mobile  charge  at  a 
io10  cur  3  concentration  near  the  surface'  is  much  more  significant  than  a 
10%  error  in  mobile  charge;  at  a  it)1**  cm  concentration  awav  from  the 
surface.  The  former  would  cause  a  nearly  10%  error  in  drain  current  while 
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the  latter  would  be  insignificant.  Actually,  the  difference  is  not  this  great. 
The  effect  is  tempered  somewhat  by  the  fact  that  as  the  charge  concentration 
decreases  away  from  the  surface,  the  gradient  decreases  also,  spreading  out 
the  lower  concentration  regions.  Thus  the  current  carried  by  the  thick  low 
concentration  regions  away  from  the  surface  is  more  nearly  equal  to  the 
current  carried  by  the  thin  high  concentration  region  near  the  surface. 

In  addition  to  inversion  layers,  accurate  representation  of  net  charge  is 
also  important  near  metallurgical  junctions.  Since  the  junctions  will  generally 
be  depleted,  the  net  impurity  concentration  is  the  driving  factor.  Although 
the  net  concentration  of  impurities  at  the  junction  is  not  high,  the  concentra¬ 
tion  gradients  near  the  junction  generally  are.  A  fine  grid  is  required,  there¬ 
fore,  in  order  to  accurately  represent  these  steep  gradients  and  to  locate  the 
junction.  Figure  2.1  shows  a  1MS0FS  grid  with  fine  spacing  normal  to  the 
surface  in  the  channel  region  and  normal  to  the  junctions  around  the  source' 
and  drain  regions. 

An  additional  point  with  regard  to  the  accurate  representation  of  charge' 
concerns  the'  allocnt  iem  eif  impurit  ies  in  a  volume  t  e>  the  imde'  representing  that 
volume.  Typically,  Urn  impurity  distribution  input  to  a  device  simulation 
program  is  evaluated  at  each  elisc re-t i zat ion  node’  ami  that  value  is  assigned 
to  the  node.  The  total  charge  within  that  mule's  volume  is  the-  proeluct 
of  the  vedume  ami  the  assignee!  impurity  concent  rat  ion.  In  regions  of  lenv 
concentration  graelicnts,  this  method  is  satisfactory;  lum-ewer,  in  re;\ie>ns  of 
high  concemtrat  iem  graelicnts,  significant  errors  in  te>tal  integrates!  charge 
may  occur.  A  better  methoel,  therfore,  is  to  integrate'  the  input  impurity 
dist  rihul  ion  over  t  he  volume  e>f  I  he'  node  t  hem  eliviele-  by  t  lie'  vedume'  te>  nrri\  v  at 
an  average  impurity  concentration.  In  this  way,  the  total  integrated  impurity 
charge-  will  he-  accurately  represent  eel.  This  method  may  allow  larger  grid 
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Fig.  2.-1.  PISCES  grid  for  an  ICFKT. 

spacing  in  the  junction  areas  of  the  device. 

2.3  Boutulnry  Sensitivities 

One  of  the  critical  steps  in  device  simulation  is  choosing  the  simulation 
region.  This  region  must  he  chosen  sullicient ly  large  so  that  the  active  region 
of  the  device  is  accurately  represented  and  is  isolated  from  the  deleterious 
effects  of  the  simulation  region  boundary.  In  general,  the  simulation  region 
must  he  large  enough  that  any  further  increase'  in  size  has  no  effect  on  the 
results  of  the  simulation.  This  is  a  check  which  may  he  used  in  practice. 
On  the  other  hand,  there  is  a  competing  desire  to  keep  the  simulation  region 
small  in  order  to  reduce  computation  time  and  memory  size.  Also,  the  hound¬ 
ary  conditions  themselves  must  accurately  represent  the  conditions  at  the 
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Boundary  condition  sensitivities  in  the  linear  region;  \\; 
V/ic;  —  Vs  —  01  ,  and  VD  =  .01 V. 
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contacts  to  the  device.  Some  device  simulation  programs  make  Amplifying 
assumptions  about  these  boundary  conditions  which  greatly  reduce  the  com¬ 
putation  time,  at  the  expense  of  reduced  accuracy  in  the  solution. 

Figure  2.5  shows  the  results  of  four  PISCFS  simulations,  each  using 
a  different  boundary  condition  equivalent  to  those  seen  in  other  simulation 
programs.  The  simulation  region  was  chosen  to  he  rectangular  in  all  four 
cases  in  order  to  maintain  equivalency.  The  device  simulated  is  an  NMOS 
FFT  with  a  metallurgical  channel  length  of  2 /mi,  gate  length  of  1/nn  (except 
Figure  2.5a)  substrate  doping  of  I0,r>  cm"3  p-type  and  an  oxide  thickness  of 

iooo  A. 

In  Figure  2.5a  the  source  and  drain  regions  are  approximated  by  vertical 
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Fig.  2.0.  Lateral  impurity  profile1  at  the  semiconductor  surface. 

boundary  contacts  1 //m  deep  as  used  in  first  generation  simulation  programs 
[2.2,  2.3] .  The  potential  applied  to  this  contact  is  reduced  slightly  from  the 
external  bias  potential  to  account  for  the  potential  drop  which  occurs  across 
the  depletion  region  in  tin1  heavily  doped  side  of  an  abrupt  junction.  In 
Figure  2.5b  the1  source  and  drain  regions  are  approximted  by  a  uniformly 
doped  rectangular  region  at  10,(1  cm  ,i.  l-’ach  region  is  1  //in  deep  and  1  //in 
long.  Figures  2.5c  and  2.5d  both  have  Gaussian  source  and  drain  profiles  with 
peak  concentrations  of  It)10  enf'  and  lateral  and  vertical  junction  depths 
of  I  //m.  The  only  difference  in  these  two  structures  is  in  the  accuracy  of 
representation  of  the  source/drain  profiles.  In  the  grid  of  Figure  2.5c.  there 
are  only  four  vertical  grid  lines  in  the  source/drain  while  in  the  grid  of  Figure 
2.5d  there  are  six.  The  continuous  doping  profile  and  the  grid  spacings  are 
shown  in  Figure  2.6. 

In  the  direct  contact  case.  Figure  2.5a,  the  oqtiipol out ial  lines  are  pulled 
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away  from  the  surface  at  the  edges.  This  means  that  the  surface  potential 
increases  near  the  source  and  drain  and  the  surface  is  more  strongly  inverted 
there  than  in  mid-channel.  Note  also  that  the  equipotential  lines  do  not 
flatten  out  in  the  channel  region  but  curve  downwards  along  the  entire 
channel  length.  This  indicates  that  the  long-channel  assumption  of  a  one- 
dimensional  potential  gradient  in  the  channel  does  not  hold  and  tin;  device 
will  exhibit  short-channel  characteristics  (i.e.  lowered  theshold  and  increased 
drain  conductance).  The  equipotential  curves  of  the  rect  angular  source/drain 
structure  are  nearly  identical. 

In  contrast,  the  Gaussian  source/drain  structure  of  Figure  2.5c  has  Hatter 
equipotential  lines  in  the  channel  region  which  do  not  pull  away  from  the 
surface  as  rapidly  at  the  channel  edges.  Since  the  surface  is  less  strongly 
inverted  the  drain  current  will  be  somewhat  less.  The  equipotential  lines 
of  the  liner  source/drain  grid  spacing  structure  of  Figure  2.5d  are  nearly 
ident  ical. 

The  drain  current  for  all  four  structures  is  plotted  as  a  function  of  gate 
voltage  in  Figure  2.7.  The  currents  were  evaluated  for  gate  voltages  of  .BY  to 
2.4V  at  ,2V  intervals  which  accounts  for  the  piecewise  continuous  appearance 
of  the  plots.  The  characteristics  of  (a)  and  (b)  are  identical  hut  differ  from 
those  of  (<•)  and  (d)  which  are  themselves  identical.  Flic  expected  variation  in 
drain  current  is  seen  to  he  approximtely  a  factor  of  two  at  large  gate  biases 
hut  nearly  an  order  of  magnit  ude  in  the  subt  heshold  region.  The  equi\ nlrnce 
of  the  drain  current  for  structures  (c)  and  (d)  implies  that  an  accurate' 
representation  of  the  source/drain  impurity  profile  may  not  he  as  important 
as  representing  it;,  general  shape,  particularly  in  the  vertical  direct  inn  near 
the  channel.  Note,  however,  in  Figure  2.t>  that  although  the  two  grid  pacings 
differ,  both  have  a  node  exactly  at  tin' junction  which  serves  to  accurately 


Fig.  2.7.  Current  sensitivity  to  boundary  conditions  in  the  subthreshohl  and 
linear  regions  for  the  structures  of  Figure  2.5. 


locate  the  junction  and  thus  fix  the  metallurgical  channel  length. 

Figure  2.8  shows  the  same  four  structures  biased  in  the  saturation  region 
of  operation.  The  pinch-off  point  is  clearly  visible  as  the  point  along  the 
channel  where  the  cquipol ential  lines  become  perpendicular  to  the  insulator- 
semiconductor  interface.  To  the  left  of  this  point,  the  electric  field  forces 
electrons  toward  the  surface  into  the  inversion  layer.  To  the  right  of  the 
pinch-off  point,  there  is  no  inversion  layer  and  the  electric  field  tends  to  spread 
out  the  electrons  as  they  travel  toward  the  drain.  In  structures  (a)  and  (b)  the 
drain  region  pushes  the  equipotentia!  lines  down  and  to  the  left  moving  the 
pinch-off  point  to  near  mid-channel  whereas  in  st  ructures  (c)  and  (d)  it  is  much 
to  the  right  of  center.  The  distance  from  the  source  to  the  pinch-off  point  is 
the  effective  channel  length  and  the  channel  length  modulation  caused  by  t  he 
pinch-off  point  moving  to  the  left  gives  rise  to  drain  conductance.  Obviously, 
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Fig.  2.8.  Boundary  condition  sensitivities  in  the  saturation  region:  I’q 
2V,  Vue  =  Vs  =  OF,  and  VD  =  5V . 


in  structures  (a)  and  (b)  the  drain  lias  greater  control  over  the  location  of  the 
pinch-off  point;  therefore,  the  currents  are  larger  and  the  drain  conductance 
is  larger  than  for  strut  ures  (r)  and  (d).  This  is  illustrated  in  Figure  2.9  where 
drain  current  is  plotted  versus  drnin-to-source  voltage  at  .5  volt  increments. 
The  gate  voltage  is  2  volts. 

In  addition  to  the  larger  currents  and  drain  conductance  expected  for 
structures  (a)  and  (l>)  it  is  also  observed  that  there  is  a  noticeable  difference 
in  current  for  (a)  and  (l>)  themselves,  (’lose  inspection  of  the  equipotential 
line  plot s  shows  that  indeed  the  pinch-off  point,  for  structure  (b)  is  slightly  to 
the  left  of  that  of  structure  (a)  resulting  in  the  larger  current.  The  source  of 
this  slight  shift,  In  vvever,  is  obvious  even  in  Fig  2.8.  If  t  In'  cquipot  ml  ini  lines 
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Fig.  2.9.  Current  sensitivity  to  boundary  conditions  in  the  saturation  region 
for  the  structures  of  Figure  2.8. 


directly  underneath  the  left  edge  of  the  drain  of  structure  (b)  arc  compared 
to  the  cquipotential  lines  underneath  the  drain  contact  of  structure  (a),  the 
oquipot.cnt.ial  lines  of  structure  (b)  are  seen  to  extend  deeper  into  the  sub¬ 
strate.  This  is  because  the  left  and  right  boundaries  oT  the  simulation  region 
of  structure  (a)  are  too  close  to  the  active  region  of  the  device  and  the  bound¬ 
ary  condition  assumptions  arc  invalid.  The  contact  portion  of  this  boundary 
is  a  fixed  potential  boundary  condition  (Dirichlet)  and  is  comparable  to  the 
rectangular  source/drain  assumption  of  structure  (b).  Along  the  remainder 
of  the  boundary  of  structure  (a),  however,  a  reflecting  boundary  condition 
(Neumann)  assumption  is  made  which  is  clearly  inaccurate.  That  is,  the 
larger  simulation  region  of  structure  (l>)  shows  that  the  cquipotential  lines 
are  not  at  all  symmetric  about  a  vertical  line  below  the  left  edge  of  the  drain. 
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This  illustrates  the  point  made  at  the  beginning  of  this  section:  the 
simulation  region  must  be  large  enough  to  be  independent  of  the  solution. 
In  particular,  reflecting  boundaries  should  only  be  placed  where  the  problem 
is  symmetric  in  some  small  region  about  the  boundary.  For  ICFETs  tlie.e 
are  three  boundaries  to  choose  since  the  top  surface  is  generally  covered  with 
the  gate  electrode  which  forms  a  natural  boundary.  The  bottom  boundary 
condition,  whether  Neumann  or  Dirichlet  must  be  placed  deeply  enough  into 
the  device  to  be  in  a  charge  neutral  region  beyond  any  depletion  regions  which 
may  arise  in  the  simulation.  Placement  of  the  left  and  right  boundaries  is 
not  as  well  defined  .  The  Dirichlet  portion  of  these  boundaries  representing 
the  source  and  drain  contact  electrodes  arc  not  critical  and  may  be  placed 
anywhere  within  the  neutral  source/drain  regions.  Th  reflecting  portions  of 
these  boundaries,  however,  dictate  that  they  be  placed  sufficiently  far  from 
the  channel  region  that  the  equipof  critin!  lines  are  r.  aurally  one  dime?:sion;>] 
(horizontal)  in  the  vicinity  of  the  boundary  for  all  bin  -  condition  -  simulate*). 
Inspection  of  Figures  2.5  and  2.8  reveals  that  none  of  the  simulation  totally 
adhere  to  this  rule  especially  along  the  drain  boundary  in  the  high  drain  bias 
case  oT  Figure  2.8.  The  simulation  region  must  contain  some  lateral  extension 
of  the  one-dimensional  portion  of  the  source  and  drain  regions  and  dace  the 
curvature  of  the  equipot enl ini  lines  extends  out  further  at  higher  drain  biases, 
the  lateral  extension  of  the  drain  region  is  significantly  larger  than  t  hat  of  the 
source  side.  Typically  one  to  two  channel  lengths  of  extension  on  the  drain 
side  is  sufficient  for  ICFETs  at.  nominal  drain  biases.  Note  that  -ince  the 
potential  and  charge  profiles  in  the  lateral  extension  regions  are  almost  one¬ 
dimensional,  a  coarse  horizontal  grid  spacing  is  sufficient  and  little  extra  grid 
is  required.  An  example  of  adequate  lateral  extension  was  shown  in  Figure 
2.1. 
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2.1  Summary 


The  various  typos  of  analysis  grids  used  for  semiconduct  or  device  simula¬ 
tion  are  presented  along  with  the  advantages  and  disadvantages  of  each. 
Methods  of  overlaying  these  grids  on  non-rectangular  structures  are  con¬ 
sidered.  The  triangular  grids  are  the  most  flexible  in  conforming  to  device 
shapes  with  the  minimum  number  of  nodes.  It  is  shown  that  the  density 
of  the  grid  should  vary  in  different  regions  of  the  device.  The  densest  grid 
should  occur  in  regions  of  high  net  charge  dnsity  or  large  gradients  of  net 
charge.  For  IGFFTs  this  means  that  the  grid  spacing  should  he  small  nor¬ 
mal  to  inversion  layers  and  metallurgical  junctions.  Checks  for  determining 
the  adequacy  of  a  particular  grid  spacing  are  suggested.  The  effect?,  of  large 
voltage  drops  between  nodes  are  presented  and  shown  to  he  highly  problem 
dependent.  An  improved  method  for  assigning  impurity  charge  to  a  node  is 
described  which  preserves  the  integrated  net  impurity  charge. 

Two  of  the  source/drain  boundary  condition  simplifications  used  br¬ 
other  device  simulation  programs  have  been  examined  and  found  to  he  grossly 
inaccurate.  Sufficient  grid  to  provide  accurate  location  of  the  source,  drain 
junctions  along  the  channel  is  necessary  but  a  coarse  lateral  spacing  in  the 
source/drain  appears  adequate  as  long  as  the  vertical  shape  is  accurately 
represented.  The  importance  of  proper  choice  of  the  simul  lion  region  has 
been  demonstrated  and  suggestions  arc  made  for  choosing;  this  space. 

The  next  chapter  describes  the  finite  difference  discio!  i.-at  ion  of  the 
model  equations  on  a  triangular  grid. 


Chapter  3 


DISCRKTIZATION 


D  sc reti/.ition  of  the  semiconductor  model  equations  using  finite- 
difTeronce  techni(|ues  on  rectangular  grids  is  a  well  established  procedure. 
Recent  work  by  Greenfield  [3.1]  has  expanded  these  procedures  to  shoe  how¬ 
to  perform  finite-difference  discretization  of  Poisson’s  equation  on  a  rectan¬ 
gular  grid  wit  h  non-planar  surface's  and  interface's.  The  use'  of  liuit ('-difference 
techniques  for  discretization  of  the'  semiconductor  mode!  equation-  on  a  tri¬ 
angular  grid,  however,  has  not  bee>n  previously  described.  The'  semiconductor 
model  equations  which  describe  the  behavior  of  semiconductor  devices  are' 
presented  in  what  follows. 

Poisson's  equation  dc'scribes  the  behavior  of  electric  llux  density  in 
regions  of  net  charge.  Since  charge's  are  sources  of  electric  llux,  the  flux 
density  must  diverge  in  regions  of  net  charge  as  given  by 


V  -D^p  1 ) 

where  I)  is  the  electric  flux  density  and  p  is  the  net  charge  concern  ration. 

The  current  continuity  equation  describes  the  time  rate  of  change  of 
carrier  concentration.  This  concentration  must  change  if  there  is  not  a 
balance  between  carrier  generation,  recombinnl  ion,  flux  into,  and  flux  out 
of  a  region  of  the  device.  'Phis  balance'  is  expressed  as 
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where  n  is  the  free  electron  concentration,  Gn  and  /?„  are  the  electron 
generation  and  recombination  rates,  q  is  the  unit  charge,  and  is  the 
electron  current  density.  The  free  hole  concentration  is  p  and  the  other 
quantities  related  to  holes  are  analogous  to  those  for  electrons. 

Carrier  transport  is  described  by  a  drift  term  dependent  upon  the  electric 
field  and  a  diffusion  term  dependent  on  the  carrier  concentration  gradient. 
This  is  given  by 


Jn  =  q[innE  +  qDnVn  (3. da) 

JP  =  quppU  -  qVp^p  (3.3b) 

where  //„  is  the  electron  mobility,  E  is  the  electric  field,  Dn  is  the  electron 
diffusion  constant  and  ;/p  and  I)p  are  the  hole  mobility  and  diffusion  constant. 

In  nondegenerate  semiconductors,  the  free  carriers  have  a  Boltzmann 
distribution  of  energy  which  leads  to  a  relationship  between  the  carrier  con¬ 
centration  and  potential  as  given  by 


(3.1a) 


p  =--  n,< 
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where  n/  is  the  intrinsic  carrier  concentration,  <p  is  the  electrostatic  potential, 
(pn  and  <pp  are  the  electron  and  hole  quasiTVrrni  potentials,  k  is  Boltzmann's 
constant  and  T  is  the  absolute  temperature.  The  nond<'g('nerat e  assumption 
is  not  required  for  this  work,  but  it  allows  simplification  and  provides  a  clearer 
picture  of  the  discretization.  Allowing  degeneracy  would  require  use  of  the 
Fermi-Dirac  distribution  function  which  leads  to  a  Fermi- Dirac  integral  for 
the  carrier  concentration  instead  of  the  exponential.  Numerical  approxima¬ 
tion  of  the  Fermi-Dirac  integral  is  computationally  no  more  expensive  than 
evaluation  of  the  exponential,  so  little  penally  is  paid  by  allowing,  degeneracy. 


Equations  (3.1)- (3.1)  represent  the  semiconductor  model.  The  following 
relations  serve  to  tie  these  equations  together.  The  net  charge  p  has  three 


components 


p  —  q{p  ~  n  +  N) 


where  N  is  the  net  ionized  impurity  concentration.  The  electrostatic  poten¬ 
tial,  electric  field,  and  electric  flux  density  are  related  by 
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where  c  is  the  permittivity.  A  link  between  the  carrier  mobility  and  diffusion 
constant  is  provided  by  the  Einstein  relation, 
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where  again  nondegeneracy  is  assumed.  Allowing  degeneracy  here  would 
require  insertion  of  a  multiplicative  factor  which  is  a  function  of  the  qunsi- 
Ecrmi  level.  Finally,  the  thermal  voltage  is  de.'iiud  as 


and  the  carrier  generation  and  reeomhinat  ion  terms  are  combined  into  a  net 
recombination  term  as 


On  =  tin  -  <’n 


vP  =  if,,  -  cP 


(3.10a) 

(3.10h) 


in  order  to  simplify  the  writing  of  equations. 

Inserting  the  relations  of  Eqs.  (3.5)-(3.10)  into  Equations  (3.1)  (3.1)  and 
assuming  steady-state  such  that  =  0,  the  semiconductor  model 

equations  become: 
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n  =  nie^'  Vr 

(3.11a) 

p  =  rue^-^l^r . 

(3.11b) 

3.1  Poisson’s  Equation 

Figure  3.1  shows  a  section  of  a  hypothetical  grid  with  five  triangular 
sections  labeled  l\  1 5  having  one  common  node  variously  labeled  t'i  i~,  and 
referred  to  as  node  i.  The  process  of  discretization  involves  the  determination 
of  two  sets  of  parameters:  the  area  assigned  to  each  node,  and  the  coupling 
coeflicients  between  pairs  of  nodes. 

The  area  assigned  to  a  node  is  taken  to  be  the  area  closer  to  that  node 
than  to  any  other  node  with  which  it  shares  a  triangle.  Tims,  in  the  live 
triangle  example  of  Figure  3.1,  the  area  /\,  hounded  by  the  da-du  d  line  I, 
represents  the  boundary  of  the  area  assigned  to  nodi'  1 .  This  boundary  is 
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[fig.  3.1.  Sample  grid  wii.li  five  triangles.  The  triangles  are  labeled  / 1  ,  .1, 

is  the  area  associated  with  the  central  node,  arid  /,  is  the  boundary 
of  that  area. 


conveniently  formed  by  the  perpendicular  bisector  of  each  edge'  common  to 
node  i. 

The  coupling  coefficient  s  are  obt  aim’d  from  t  he  disc  ret  i/a  l  ion  of  1’oisaon’s 
equation.  This  is  achieved  by  applying;  Gauss’  law  to  Ifip  (3.1)  ;  ad  converting 
integrations  into  summations.  Applying;  Gauss’  law  to  lap  (3.1), 


I)  ■  dl 


Inserting  the  relation  of  lap  (3.7)  and  recognizing  that  (he  integrals  may  be 
evaluated  by  parts  common  to  each  triangle, 
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where  /im  and  Airn  are  those  portions  of  the  node  t  boundary  and  area  lying 
within  triangle  tm. 

At  this  point,  two  assumptions  are  made:  the  permittivity  and  the 
electric  field  are  constant  within  a  triangle.  The  first  of  these  requires  only 
that  material  boundaries  lie  along  triangle  edges,  a  requirement  easily  met 
with  the  flexible  triangular  grids.  The  assumption  of  constant  electric  field 
within  a  triangle  is  a  natural  result  of  the  triangular  discretization  since  the 
only  unique  interpolation  function  of  three  points  in  two  dimensions  is  a  linear 
function,  and  linear  potential  implies  a  constant  electric  field. 

These  assumptions  allow  the  line  integrals  to  be  replaced  by  dot  products. 
Furthermore,  since  the  net  charge  assigned  to  a  node  is  considered  to  be 
evenly  distributed  throughout  its  area,  p  is  a  spacial  constant  and  may  be 
taken  out  of  the  integral.  These  changes  to  Eq.  (3.16)  yield 

‘  (jiljl  +  ^ilk\ )  +  f-lE  •  (ji2j2  +  •  •  ’  +  C5 IC  ■  (jibjj  +  T, r, /c 5 ^ 

=  PA-  (3.17' 

where  /tmjm  is  the  vector  normal  to  the  portion  of  the  boundary  /,  in  triangle 
tm  which  is  perpendicular  to  the  ij  side.  The  magnitude  of  is  equal  to 

the  length  of  the  boundary  segment  and  the  positive  direction  is  away  from 
node  z.  The  other  boundary  normal  vectors  are  similarly  defined. 

Figure  3.2  shows  the  labeling  convention  for  a  sample  triangle.  The 
vectors  Uj  and  are  the  unit  vectors  in  the  ij  and  ik  directions  respectively. 
The  subscript  denoting  the  triangle  number  has  been  dropped  in  the  figure 
and  the  boundary  segments  lij  and  have  been  relabeled  with  their  length 
hj  and  h k.  Equation  (3.17)  may  then  be  written  in  summation  notation  as 

^  ^  '  (hj’m^jm  T  ^-km^km)  =  (3.  IN) 
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Fig.  3.2.  Labeling  convention  for  a  triangle.  The  distance  between  nodes  is 
d,  the  length  (height)  of  the  boundary  segments  is  h,  and  u  is  a 
unit  vector. 


where  the  problem  has  been  generalized  to  include  any  number  A/  oT triangles 
containing  the  common  node  i. 

The  dot  product  E  •  Uj  is  th<;  component  oF  the  electric  field  in  the  ij 
direction.  Due  to  the  assumption  of  constant  electric  field  within  t  he  triangle, 
this  may  be  discretized  using  Equation  (3.G)  to  obtain 


E  •  Uj  =  — 


dj 


(3.19) 


where  dj  is  the  distance  between  nodes  i  and  j.  Equation  (3.18)  may  then  be 
written  as 


Yl  Crn((^x  -  m)h,:rn  +  [4>i  ~  =  PiA-  (3-20) 
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This  is  the  discretized  form  of  Poisson’s  equation.  Summing  the  terms  over 
all  triangles  containing  node  i  results  in  one  equation  containing  the  unknown 
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potentials  of  node  i  and  all  adjoining  nodes.  The  term  js  referred  to 

as  the  coupling  coefficient  between  nodes  i  and  j  in  triangle  m.  There  is  a 
similar  coefficient  for  the  same  two  nodes  in  the  adjacent  triangle  (if  there  is 
one)  which  adds  to  this  to  form  the  full  coupling  coefficient.  Note  in  liquation 
(3.20)  that  the  coefficient  of  every  node  adjacent  to  node  i  is  also  a  coefficient 
of  node  i  with  opposite  sign,  thus  the  coefficient  of  ipi  is  positive  (r,  h  and  d 
are  all  positive)  and  exactly  equal  to  the  negative  of  the  sum  of  the  coefficients 
of  all  of  the  adjacent  nodes. 

Repeating  the  summation  of  Eq.  (3.20)  for  each  node  in  the  grid  results 
in  a  set  of  N  equations  in  N  unknowns  where  N  is  the  number  of  nodes  in  t  he 
grid.  The  coefficient  matrix  for  this  set  of  equations  is  diagonally  dominant. 
Note  that  if  the  grid  is  rectangular  (properly  composed  of  right  triangles)  the 
set  of  equations  reduces  to  exactly  that  of  the  standard  five-point  difference 
scheme. 

Unfortunately,  the  charge  concentration  pi  is  a  function  of  potential  as 
described  in  Eq.  (3.14).  Combining  Eqs.  (3.5)  and  (3.1  1)  for  node  i  gives 


q(n,el+*-*'VVT  -  nIct+‘-+'-'VvT 


(3.21) 


thus  the  potential  V\  appears  non  linearly  in  the  right  hand  side  of  Eq.  (3.20). 
The  resulting  non  linear  equation  is  solved  using  Newton’s  method.  That  is, 
it  is  linearized  and  solved  iteratively  until  converged. 

Combining  equations  (3.20)  and  (3.21),  linearization  is  .achieved  by  first 
replacing  every  V;  by  +  A»/’, 
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(3.22) 


Y1  fm(  ((i>i  +  A^,-)  -  {ipjm  +  A 

l<m<M  '  d>m 

+  ((V’i  +  AV»i)  -  (V’fcm  + 

«fcm  / 

=  qAi(nlel*'i"l*i  +  ^*i))/VT  -  n/e((v'1  +  A’/'.)-<^n,)/ vT  +  yy. 

=  <Ml{p.e_AV’*/V'7'  -  n,e^^VT  +  JV.) 
then  taking  first  order  Taylor-series  approximations  to  the  exponentials, 

=  q/l,(pi(l  -  A ilu/Vr)  ~  n,(l  +  A 4’i/Vr)  +  1 Vt) 

=  P,A>  -  9.Ai(p,-  +  n^Ai/’i/Vr.  (3.23) 


Placing  all  terms  in  At/’  on  the  left  hand  side  results  in  the  iterative  form, 


y;  (  m  ( (  A  )/>,  -  A  l !'j  m )  hJ  m  +  ( A  0 ,•  -  A  l/’fc m )  y 

V  dJ - 


1  <  rn  <  A/ 


+  7-A1(p,  +  n^&ipjVr 


--  E 


1  <  m  <  M 


(01  “  0jm)  +  (0j  -  0fcm)  y”‘  I  +  /V'i-  (3.21) 


The  right  hand  side  of  Fq.  (3. 21)  is  the  residual  of  Poisson’s  equation  and  (lie 
left  hand  side  contains  the  unknown  potentials  for  the  Newton  step.  After 
each  solution  for  the  A V’’s,  the  potentials  are  updated,  the  potential  depen¬ 
dent  terms  are  re-evaluated,  and  the  iteration  proceeds  until  convergence  is 
achieved. 

When  Fq.  (3.2  1)  is  assembled  for  every  point  in  the  grid  and  put  in  matrix 
form,  the  coefficient  matrix  has  exactly  the  same  terms  in  it  as  the  coefficient 
mat  rix  of  Fq.  (3.20)  except,  that  positive'  terms  resulting  from  the  linearization 
of  the  charge  concentration  have  been  added  to  the  entries  on  the  diagonal. 
This  matrix  is  positive  definite  resulting  in  desirable  convergence'  properties. 


3.1.1  Area  Allocation 


In  the  foregoing  discussion,  the  vectors  hJ„luJTn  and  of  Kq.  (3.18) 

were  kept  separate  in  order  that  the  electric  field  may  be  expressed  in  terms 
of  its  components  in  the  ujni  and  directions.  As  shown  in  ICq.  (3.20) 

these  components  are  conveniently  obtained  from  the  potential  differences 
between  nodes.  Additional  insight  may  be  gained,  however,  by  performing 
the  vector  sum  as  illustrated  in  f  igure  3.3.  The  electric  flux  crossing  the  node 
i  area  boundary  l,  within  this  triangle  is  identically  the  Ilux  crossing;  the  line 
segment  s  joining  the  midpoints  of  the  adjacent  sides,  thus  the  choice  of  the 
node  i  area  boundary  is  somewhat  arbitrary.  In  fact,  the  above  discretization 
applies  to  any  simple  boundary  which  has  these  midpoints  as  its  endpoints. 
Four  random  possibilities  are  shown  in  Figure  3.1. 

One  such  alternative  boundary  is  to  use  the  lino  segments  joining'  the 
triangle  centroid  to  the  adjacent  side  midpoints  in  the  manner  of  Winslow 
[3.2].  With  this  boundary,  each  node  is  allocated  one-third  oT  the  total 
triangle  area.  This  allocation  scheme  was  tried  with  the  F1SCFS  program  and 
generally  yielded  satisfactory  results  except  that  some  assemblies  oT  triangles 
resulted  in  uneven  distribution  of  area.  An  example  is  shown  in  Figure  3.5 
in  which  two  similar  discretization  grids  arc  given  with  the  only  difference 
being  the  triangle  orientation.  Note  that  in  case  (a)  there  are  four  identical 
triangles  connected  to  the  central  node  while  in  case  (b)  there  are  eight. 
As  a  result,  the  central  node  is  allocated  twice  as  much  area  in  case  (b)  as 
in  case  (a).  This  unequal  weighting  results  in  perturbations  to  the  desired 
solution.  In  contrast,  application  of  the  perpendicular  bisector  method  to  the 
grids  of  Figure  3.5  results  in  equal  allocation  of  area  for  the  two  cases.  The 
perpendicular  bisector  method  of  forming  the  node  area  boundary  appears 
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Fig.  3.5.  Dependence  of  area  weighting  on  triangle  orientation  using  the 
centroid  method.  The  area  allocation  in  case  (b)  is  twice  that  of 
case  (a). 

to  be  the,  best  area  allocation  scheme  for  triangular  grids  in  semiconductor 
problems. 

3.1.2  Obtuse  Triangles 

The  above  derivations  work  quite  nicely  for  acute  triangles;  however, 
when  the  triangles  become  obtuse,  adjustmets  need  to  be  made.  The  need 
arises  from  the  fact  that  tiie  intersection  of  the  perpendicular  bisectors  of 
the  sides  of  an  obtuse  triangle  occurs  outside  of  the  triangle  as  illustrated  in 
Figure  3.G.  In  this  example,  the  coupling  coefficients  are  exactly  as  derived 
earlier  except  that  the  coupling  coediciont  between  nodes  i  and  j  becomes 
negative,  —  since  the  segment  hj  has  reversed  direction.  It  can  be  shown 
that  this  choice  of  coupling  coefficients  still  accurately  accounts  for  the  flux 
passing  through  the  line  segment  *.  Intuitively,  it  may  be  reasoned  that  the 
vector  sum  of  the  directed  line  segment  /ifc  minus  the  directed  line  segment  h j 
is  the  directed  line  segment  s,  so  the  vector  sum  of  their  normals  also  equate. 
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Fig.  3.6.  Coupling  coefficients  for  Poisson’s  equation  on  an  obtuse  triangle. 


Note  that  the  coupling  coefficient  goes  from  positive  to  zero  to  negative  as 
the  triangle  goes  from  acute  to  right  to  obtuse,  so  there  is  no  discontinuity 
in  the  transition  from  acute  to  obtuse. 

The  total  coupling  coefficient  for  nodes  i  and  j  includes  a  contribution 
from  both  triangles  which  share  the  side  joining  nodes  i  and  j.  Thus,  ah  hough 
one  component  of  the  coupling  coefficient  may  be  negative,  the  total  coupling 
coefficient  may  still  he  positive  if  the  contribution  from  the  adjacent  triangle 
is  sufficiently  large  and  positive.  It  can  be  shown  that  the  total  coupling 
coefficient,  will  be  positive  if  the  sum  of  the  opposite  angles  is  less  than  1*0 
degrees.  This  condition  can  always  be  satisfied  for  triangles  which  are  not  on 
the  solution  region  boundary  or  on  a  material  interface  by  reconnecting  the 
triangles.  That  this  is  so  can  easily  be  proven  by  hypothesizing  a  situation  in 
which  two  adjacent  triangles  exist  whose  angles  opposite  the  common  edge 
sum  to  more  than  180  degrees.  If  the  four  nodes  composing  these  two  triangles 
are  looked  upon  as  a  quadrilateral,  then  the  common  edge  is  a  diagonal  of 
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the  quadrilateral,  the  sum  of  whose  interior  angles  must  equal  360  degrees. 
If  instead  of  dividing  the  quadrilateral  into  two  triangles  using  this  diagonal, 
the  other  diagonal  is  used,  two  triangles  result  whose  angles  opposite  their 
common  side  necessarily  sum  to  less  than  180  degrees.  Obviously,  if  the  edge 
in  question  composes  a  solution  region  boundary  or  a  material  interface  ,  the 
triangles  may  not  be  reconstructed;  however,  subdivision  of  a  triangle  into 
two  or  more  triangles  can  eliminate  the  problem  in  these  special  cases. 

The  occurrence  of  obtuse  angles  also  complicates  the  area  allocation 
method.  Looking  again  at  Figure  3.6,  the  boundary  for  the  area  allocated  to 
node  i  is  no  longer  defined  by  the  line  segments  hj  and  h/<  as  they  were  in 
Figure  3.2.  One  possible  choice  of  boundary  is  that  portion  of  line  segment 
which  lies  within  the  triangle.  Unfortunately,  this  choice  does  not  meet  the 
prerequisite  that  the  boundary  include  the  midpoints  of  the  adjacent  sides, 
thus  there  would  not  be  conservation  of  flux  within  the  triangle  using  this 
boundary. 

A  better  choice  for  the  boundary  of  the  area  allocated  to  node  i  is 
the  line  segment  s  in  Figure  3.6,  the  segment  joining  the  midpoints  of  the 
adjacent,  sides.  This  choice  satisfies  the  condition  for  conservation  of  flux 
and  is  somewhat  better  than  the  centroid  method  in  area  weighting  in  that 
it  allocates  one-fourth  of  the  triangle  area  to  node  one-fourth  to  node  j, 
and  one-half  to  node  k.  The  acute  and  obtuse  area  allocation  schemes  are 
identical  for  an  angle  of  90  degrees  so  there  is  a  smooth  transition  from  one 
to  the  other.  Figure  3.7  shows  a  hypothetical  grid  with  various  triangle  types 
to  demonstrate  the  area  allocation  scheme.  The  solid  lines  are  the  triangular 
grid  and  the  dashed  lines  are  the  area  boundaries. 


Fig.  3.7.  Area  allocation  with  a  mix  of  acute  and  obtuse  triangles. 


3.2  Continuity  Equation 

The  discretization  of  the  continuity  equation  will  be  described  only  for 
electrons  since  the  discretization  for  holes  is  analogous  with  the  exception 
of  signs.  Also,  the  discretization  of  the  electron  transport  equation  will  be 
performed  in  one  dimension  only.  The  reason  for  this  will  be  explained  later. 


3.2.1  Electron  Transport  Equation 


From  Eq.  (3.13),  the  electron  transport  equation  in  one  dimension  is 

d 


Jn  —  <7/tn(  nE  T  V'l 


(3.25) 


The  standard  finite-difference  approach  to  discretization  of  this  equation 
would  result  in 


./»( V)  =  9/‘n(V)^n(  V)^(‘V)  +  VV 


rt{A  .r)  —  n(0) 


Ax 


) 


(3.2B) 
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where  a  uniform  spa  ing  of  Ax  is  assumed  for  simplicity.  Insertion  of  Eq. 
(3.28)  into  the  continuity  equation  (3.12a)  yields  a  single  equation  in  the  un¬ 
known  n  assuming  that  the  electric  field  is  given.  However,  it  can  be  shown 
that  this  form  of  the  transport  equation  leads  to  numerical  instability  when¬ 
ever  the  voltage  change  between  nodes  exceeds  2VV .  This  point  was  noted  by 
Scharfettcr  and  Gummel  [3.3]  and  has  resulted  in  formulation  of  a  discretiza¬ 
tion  scheme  which  avoids  this  difficulty.  The  method  will  he  summarized 
here  in  order  to  provide  needed  background  for  further  discussion. 

Equation  (3.25)  may  be  rearranged  as 


dn  E  Jn 

- f- - n  =r - — - 

dx  VT  QUnVr 


(3.27) 


Given  two  nodes  separated  by  a  distance  d  as  in  Figure  3.8  and  assuming 
that  Jn,  nn>  and  E  are  constant  between  these  nodes,  then  Eq.  (3.27)  may 
be  viewed  as  a  first  order  differential  equation  in  n  with  constant  coefficients. 
The  general  solution  to  this  equation  is 


n  =  C'-B’lv*  + 

qHnE 

Evaluating  this  equation  for  n  ~  n,  at  x  =  0  results  in 

Jr, 


C  =  Tli 


qHnE 


Using  this  value  for  C  and  evaluating  Eq.  (3.28)  at  node  j  yields 

Jr, 

qn„E 


V-'w/’v  + 

'  qfinEJ 


which  may  be  rearranged  as 


r  rrf  7lJ  ni  \ 

n  —  <7/'n  _e- /■;«// vT  +  j  _cm/vT  ) 


(3.28) 


(3.29) 


(3.30) 


(3.31) 
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Fig.  3.8.  Labeling  convention  for  one-dimensional  current  transport. 


or  equivalently 


r  _ 

Jn~Wn  d 


+ 


n, _  \ 

1  _  g(V>t  —  Vo)/ w  ) 


(3.32) 


This  is  the  Seharfetter-Gumme!  form  for  discretization  of  the  electron  trans¬ 
port  equation  in  one  dimension.  One  may  easily  verify  that  in  the  limit  as 
the  potential  difference  between  nodes  approaches  zero,  Fq.  3.32  becomes 
identically  the  diffusion  term  while  as  the  potential  difference  becomes  large, 
it  becomes  identically  the  drift  term. 


3.2.2  Electron  Continuity  Equation 

The  discretization  of  the  electron  continuity  equation  will  be  performed 
in  a  quasi- two-dimensional  form  using  the  one-dimensional  t  ransport  equa¬ 
tion.  As  with  Poisson’s  equation,  we  begin  with  the  acute  triangle  case  and 
use  the  boundaries  shown  in  Figures  3.1  and  3.2.  Applying  Gauss’  Theorem 
to  Eq.  (3.12a), 


qUndA. 


(3.33) 
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The  discretization  process  from  this  point  parallels  that  of  the  Poisson  equa¬ 
tion  (Eq.  3.15  to  Eq.  3.20)  except  that  quantities  on  the  left  hand  side  of  the 
equation  are  taken  to  be  constant  only  along  an  edge  of  the  triangle  instead 
of  over  the  entire  triangle.  The  discretization  leads  to  a  form  which  is  similar 
to  that  of  Eq.  (3.18) 

^  ^  ^Jnjm  "  d~  Jnkm  *  q^J-n^i  (3.34) 

1  <  m  <  M  V 

where  Jnjrn  is  the  electron  current  density  along  edge  ij  in  triangle  m  and 
Jnkm  is  similarly  defined.  Note  that  since  Jn  cannot  be  assumed  constant 
over  the  triangle,  Jnjm  7^  Jnkm ■  Since  Jnjrn  and  Jnk,n  a™  one-dimensional 
current  density  vectors  in  the  Ujm  and  directions  respectively,  the  dot 
products  become  multiplies  and  Eq.  3.34  becomes 

y  ^  (Jnjmhjm  d"  Jnkm^km)  Q  /l ,  (3.3o) 

l<m<  M 

where  the  current  density  is  taken  to  be  positive  in  the  direction  away  from 
node  i.  Thus  the  electron  current  out  of  the  area  of  node  i  into  the  area  of 
node  j  in  triangle  m  is  the  current  density,  Jnjm ,  times  the  length  of  the 
perpendicular  boundary  through  which  the  current  Slows,  hjm.  This  same 
current  density  also  flows  through  a  perpendicular  segment  of  the  node  i 
boundary  in  the  adjacent  triangle.  The  total  electron  current  out  of  the  node 
i  area  is  the  summation  of  the  current  crossing  the  area  boundary  into  each 
adjacent  node,  and  this  total  must  equal  the  recombination  rate  qUn.\x. 

Finally,  Eq.  (3.32)  may  be  inserted  into  Eq.  (3.35)  dropping  the  q  from 
both  sides  to  yield 
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l<m<MV  Jm  e 


+ 


n. 


)/vT  1  _  e(ri>i-*jm)/VT 


) 


,  ^nkm^km  /  ,  ,  \ 

+  - ,  - - (V\  -  V»Jfcm) 

d-km 

=  C/„Ai. 


^ _ nkm _  ^  t _ \  \ 

1  —  e~^'~^krn)/VT  1  —  e^'~^kn‘^vT  ) ) 


(3.36) 


This  is  the  fully  discretized  electron  continuity  equation.  The  unknowns  are 
the  electron  concentrations  n,,  njrn,  and  rifcm  while  the  mobilities,  potentials, 
and  electron  recombination  rate  are  assumed  to  be  known;  although,  each 
is  actually  a  function  of  the  electron  concentration.  The  potentials  are  a 
strong  function  of  the  electron  concentration  through  Poisson’s  equation,  thus 
alternating  solutions  of  the  Poisson  and  continuity  equations  are  required 
until  convergence  is  obtained.  The  mobilites  and  electron  recombination 
rates  may  also  be  functions  of  the  electron  concentrations,  but  the  functonal 
relationships  are  generally  very  weak  so  that  merely  updating  the  mobility 
and  recombination  rate  after  new  electron  concentrations  arc  computed  is 
sufficient. 

When  Eq.  (3.36)  is  assembled  for  every  node  in  the  grid,  one  again  has 
N  equations  in  N  unknowns  as  in  the  Poisson  case;  however,  the  continuity 
equation  docs  not  have  to  be  solved  in  the  insulator  regions  of  the  device  so 
that  one  may  limit  the  solution  to  only  those  nodes  lying  in  the  semiconductor 
region.  The  coefficient  matrix  does  not  have  the  desirable  iteration  properties 
of  the  Poisson  coefficient  matrix,  but  since  iteration  is  not  required  to  solve 
the  continuity  equation,  this  point  is  not  critical.  Also,  the  coefficients  look 
rather  formidable  at  first,  but  a  careful  look  will  show  that  the  coefficients 
are  well  behaved  and,  in  fact,  no  difficulties  in  obtaining  a  solution  were  ever 
observed  in  the  PISCES  program.  The  well  behaved  nature  of  the  coefficients 
is  demonstrated  in  the  fact  that 
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(3.37) 


Finally,  note  that  the  coupling  terms  in  Eq.  (3.36)  are  analogous  to  the 
terms  'j1  in  Poisson’s  equation,  Eq.  (3.20). 


3.2.3  Obtuse  Triangles 

As  in  the  case  of  Poisson’s  equation,  adjustments  need  to  bo  made  to 
the  continuity  equation  discretization  in  obtuse  triangles.  In  the  Poisson 
dciscrctization,  the  node  area  boundary  was  pinned  at  the  midpoints  of  the 
sides  as  the  triangle  became  obtuse  while  the  coupling  coefficients,  -■} ,  were 
allowed  to  become  negative  in  a  rather  natural  manner.  The  continuity 
equation  necessarily  uses  the  same  boundary  as  Poisson’s  equation  but  cannot 
allow  its  -j-  terms  to  become  negative  as  the  trangle  becomes  obtuse.  This 
discrepancy  stems  from  the  quasi-two-dimensonal  nature  of  the  continuity 
equation  discretization. 

Figure  3.9  graphically  depicts  the  reasoning  behind  the  choice  of  coupling 
coefficients  for  the  continuity  equation  in  acute  and  obtuse  triangles.  The 
current  density  flowing  from  node  a  node  b  is  shown  as  Jn.  In  triangle  l-i  (bis 
current  density  passes  from  the  node  j'a  area,  into  the  node  area  through 
a  cross-sectional  window  of  width  hj2  which  is  the  length  of  the  boundary 
segment.  (Note  that  node  a  is  equivalent  to  both  in  tl  and  i%  in  1 2). 
However,  in  the  obtuse  triangle,  t  j,  the  current  density  passes  through  a 
cross-sectional  width  of  h^i  which  is  the  projected  length  of  the  boundary 
segment  onto  the  perpendicular  to  the  current  density  vector.  The  value  of 


55 


Fig.  3.9.  Coupling  coefficient  derivation  for  the  continuity  equation  on  acute 
and  obtuse  triangles. 

h  used  in  the  total  coupling  coefficient  between  nodes  a  and  b  is  the  sum  of 
these  two  terms,  i.e.  the  total  cross-sectional  width  of  the  common  boundary 
between  the  two  nodes. 

Note  in  obtuse  triangle  t\,  that  there  is  no  common  boundary  between 
nodes  i\  and  j\  so  hji  —  0  and  no  current  is  allowed  to  flow  between  them. 
Current  may  be  allowed  along  this  same  edge  in  an  adjacent  triangle,  however, 
if  the  opposite  angle  is  acute.  This  condition  can  be  assured,  as  was  shown  in 
Section  3.1.2,  as  long  as  the  triangle  edge  is  not  a  solution  region  boundary 
or  a  material  interface.  Even  these  cases  can  be  handled  by  subdivision  of 
the  triangles,  so  it  is  never  necessary  for  the  continuity  equation  coupling 
coefficient  between  two  adjacent  semiconductor  nodes  to  be  zero  (except  in  the 
degenerate  case  of  a  perfectly  rectangular  grid  structure  where  the  triangle 
edges  representing  the  diagonals  of  the  rectangles  have  coefficients  of  zero). 
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Note  also  that  as  in  the  Poisson  discretization,  there  is  no  discontinuity  in 
the  choice  of  coupling  coefficients  as  the  triangle  passes  from  acute  to  obtuse. 


3.2.4  Absence  of  Two-Dimensional  Scharfetter-Gummel  Discretization 


As  mentioned  earlier,  the  Scharfetter-Gummel  algorithm  for  discretiza¬ 
tion  of  the  electron  transport  equation  does  not  have  an  analogous  form  in 
two  dimensions.  This  can  be  shown  as  follows.  Equation  (3.13a)  may  be 
rewritten  in  two  dimensional  form  as 

J n  4/^ u v x  -{  Ayiiy)  -f-  V 7’ ^  ^ ^  ilx  +  Qy  ^ y ^ ^ 
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where  ux  and  uy  are  the  unit  vectors  in  the  x  and  y  directions.  Making  the 
necessary  assumptions  that  /t„,  Jn  and  E  are  constant  over  the  triangle,  then 
the  x  and  y  components  of  lap  (3.38)  may  be  handled  separately  using  the 
one-dimensional  algorithm.  In  the  x  direction, 

Bix  /tj  %j ux 

Bx  VTn  +  <iy,  ,VT 

—  mi  -\-  p  (3.30) 

where  n  and  fi  take  on  the  obvious  values.  Similarly,  in  the  y  direction, 

Bxi  I[/y  B ny 
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—  7  ti  +  <5.  (3.10) 

It  is  a  necessary  and  sufficient  condition  for  a  solution  to  exist  for  Eqs. 
(3.30)  and  (3.40)  that 

<)2n  d2n 
()x<hj  dy  dx 
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Evaluating  these  terms  yields 


and 


d2n 

dxdy 


=  Tx[ln  +  s) 

=-----  7( an  4-  /3) 


d2n 

dydx 


=  ~(an  +  fi) 
—  a^n  +  6). 


(3.42) 


(3.43) 


Thus,  the  necessary  and  sufficient  condition  for  a  solution  to  exist  is 


or 
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Inserting  the  values  of  a,  /?,  7,  and  6  from  Eqs.  (3.39)  and  (3.10), 
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Therefore,  a  solution  exists  if  and  only  if  E  and  J  are  collincar;  that  is,  only 
if  the  problem  is  one-dimensional. 


3.3  Summary 

The  model  equations  were  introduced  and  finite-difference  discretization 
on  an  irregular  triangular  grid  was  described.  Discretization  of  Poisson’s 
equation  on  an  acute  triangle  was  shown  and  subtleties  of  area  weighting 


58 


schemes  were  addressed  with  the  perpendicular  bisector  method  shown  to  be 
preferable.  Modifications  to  the  area  weighting  were  shown  to  be  necessary 
for  obtuse  triangles.  Also,  coupling  coefficients  within  individual  triangles 
were  seen  to  be  negative  for  the  edge  opposite  an  obtuse  angle  to  achieve 
proper  conservation  of  electric  flux.  However,  it  was  demonstrated  that  the 
total  coupling  coefficient  (the  sum  from  two  adjacent  triangles)  need  never 
be  negative  as  long  as  the  triangle  edge  opposite  the  obtuse  angle  is  not  a 
solution  region  boundary  or  an  interface  between  dissimilar  materials. 

Discretization  of  the  electron  continuity  equation  was  described  in  a 
quasi-two-diincnsonal  form  using  the  Scharfctter-Gummel  algorithm  for  dis¬ 
cretization  of  the  electron  transport  equation.  Obtuse  triangles  were  also 
shown  to  require  special  consideration  in  computing  the  coupling  coefficients 
for  the  continuity  equation  discretization.  The  absence  of  a  two-dimensional 
form  of  the  Scharfetter-Gurnmel  algorithm  was  demonstrated. 

The  next  chapter  covers  the  various  methods  of  solving  the  combined 
Poisson  and  continuity  discretized  equations,  factors  affecting  convergence, 
and  means  of  accelerating  the  convergence. 
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Chapter  4 


SOLUTION  TECHNIQUES 

Discretization  of  Poisson’s  equation  and  the  continuity  equation  results 
in  two  systems  of  equations  which  must  be  solved  in  order  to  determine 
the  potentials  and  carrier  concentrations.  There  are  two  aspects  to  the 
solution  of  the  equations— one  is  the  solution  of  the  matrix  equation  A/x  = 
y  representing  either  the  discretized  Poisson  or  continuity  equation  by  itself, 
the  other  is  the  consistent  solution  of  the  two  coupled  equations  together. 
This  chapter  deals  with  both  of  these  aspects. 


4.1  Matrix  Equation  Solution  Methods 

In  this  section,  various  grid  solution  method  combinations  are  discussed 
for  both  Poisson’s  equation  and  the  continuity  equation.  It  is  important  to 
note  that  the  coefficient  matrices  in  these  matrix  equations  are  sparse  [4.1] 
that  is,  most  of  the  matrix  elements  are  zero.  The  matrix  elements  are 
the  coupling  coefficients  between  nodes  in  the  grid,  so  only  those  elements 
corresponding  to  coupled  nodes  will  be  non-zero.  This  sparsity  must  be 
exploited  to  minimize  the  equation  solution  time. 

The  matrix  equation  solution  methods  may  be  divided  into  two  classes: 
direct  and  iterative.  In  the  direct  solution  techniques,  a  linear  equation  solu¬ 
tion  is  achieved  in  a  deterministic  number  of  steps  of  the  solution  algorit  hm. 
The  number  of  steps  depends  only  on  the  algorithm  chosen  and  the  connec¬ 
tivity  of  the  grid  (i.c.  how  the  nodes  are  interconnected)  and  is  independent 
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of  the  values  of  the  coupling  coefficients  or  solution.  Iterative  matrix  solu¬ 
tion  techniques  have  a  deterministic  number  of  steps  for  one  pass  (iteration) 
through  the  algorithm,  but  several  or  many  iterations  of  the  algorithm  must 
be  made  in  order  to  obtain  an  accurate  solution,  thus  the  total  number  of 
steps  required  is  non-dcterministic.  In  particular,  the  number  of  iterations 
depends  on  the  values  of  the  coupling  coefficients,  the  values  of  the  solution, 
and  the  accuracy  of  solution  desired. 

Typically,  the  number  of  operations  (multiplies  and  adds)  for  a  direct 
solution  is  significantly  larger  than  for  a  single  pass  through  an  iterative 
solution.  As  the  number  of  iterations  increases,  however,  the  tot  al  number  of 
operations  required  for  a  solution  using  an  iterative  algorithm  may  surpass 
that  required  for  a  direct  solution,  thus  there  is  a  break-even  point  where  one 
method  becomes  more  efficient  than  the  other.  Unfortunately,  locating  the 

' 

break-even  point  is  more  of  an  art  than  a  science, 

f 

4.1.1  Poisson’s  Equation 


As  mentioned  in  Chapter  3,  the  coefficient  matrix  for  Poisson’s  equation 
is  symmetric  am,  positive  definite.  The  symmetry  may  be  exploited  to  reduce 
coefficient  storage  and  operation  count.  The  positive  definite  characteristic 
is  required  for  convergence  of  some  of  the  iterative  techniques. 

Of  the  direct  solution  methods,  LU  decomposition  [1.2]  is  the  most 
common.  Based  on  Gaussian  elimination,  it  may  be  used  to  solve  full  or 
sparse  matrix  equations.  Sparse  LU  decomposition  merely  ignores  the  zero 
entries  in  the  coefficient  matrix  and,  in  this  sense,  may  be  considered  a  'brute 
force”  method.  The  coefficient  matrix,  A / ,  is  decomposed  into  the  product, 
of  a  lower  triangular  matrix  L  and  an  upper  triangular  matrix  l 1  such  that 


61 


M  —  LU .  The  matrix  equation  LUx  =  y  is  then  easily  solved  by  using 
forward  substitution  on  Lw  =  y  and  backward  substitution  on  Ux  =  w. 

The  advantages  of  this  method  are  its  flexibility  and  stability.  It  can 
be  used  on  literally  any  grid,  with  any  discretization,  and  with  almost  any 
non-singular  coefficient  values.  Since  the  PISCES  program  was  written  in 
order  to  study  grid,  discretization,  various  device  types,  and  various  regions 
of  operation,  this  flexibility  and  stability  was  extremely  desireable.  The 
disadvantages  of  LU  decomposition  are  that  it  requires  large  amounts  of 
storage  and  has  a  high  operation  count  for  large  grids.  The  large  storage 
arises  from  the  fact  that  the  L  and  U  matrices  both  have  more  non-zero  entries 
than  the  coefficient  matrix  A/,  a  condition  known  as  fill.  This  also  influences 
the  operation  count  which  is  on  the  order  of  N2 ,  (written  as  0(.Y“)),  where 
N  is  the  number  of  nodes  in  the  grid. 

The  second  direct  solution  method  of  interest  is  the  Fast  Fourier 
Transform  (FFT)  method  as  applied  to  the  solution  of  Poisson’s  equation  by 
Hockney  ['1.3]  for  plasma  physics  computations.  It  is  quite  the  opposite  of  LU 
decomposition  in  the  sense  that  it  has  small  storage  requirements  and  a  low 
operation  count,  but  is  very  inflexible  in  terms  of  grid  or  device  type.  The 
FFT  method  demands  grids  similar  to  those  of  Figures  2.01b  or  2. Old  that 
is,  the  horizontal  grid  spacing  must  be  uniform  along  any  given  line  and  the 
number  of  nodes  along  that  line  should  be  a  power  of  two.  The  method  is 
based  on  the  orthogonality  of  Fourier  harmonics  -  if  a  function  is  an  electros¬ 
tatic  potential  solution  of  Poisson’s  equation  for  a  given  charge  distribution, 
then  each  Fourier  harmonic  of  the  function  is  a  solution  of  Poisson’s  equa¬ 
tion  for  the  corresponding  harmonic  of  the  charge  distribution.  The  method 
works  most  efficiently  if  the  Fourier  analysis  is  performed  in  one  dimension 
only  which  is  normally  the  lateral  direction  for  semiconductor  devices. 
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In  the  discrete  problem  there  are  only  as  many  harmonics  as  there  are  nodes 
in  the  x  direction  along  the  grid  line,  so  for  an  m  by  n  grid  the  problem 
has  been  reduced  from  a  two-dimensional  matrix  equation  of  oider  N  (where 
Ar  =  m?i)  to  n  one-dimensional  matrix  equations  of  order  rn.  Since  the 
matrix  equations  are  one-dimensional,  the  coefficient  matrix  is  tridiagonal 
and  very  efficient  solution  methods  can  be  used.  The  total  operation  count 
is  0(  iV  In  rn). 

There  are  several  disadvantages  to  the  use  of  this  method.  One  is  the 
strict  requirement  for  a  rectangular  grid  with  uniform  horizontal  grid  spac¬ 
ing.  Another  is  the  requirement  that  permittivity  be  constant  in  the  lateral 
direction,  limiting  its  application  to  devices  with  planar  material  interfaces 
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and  surfaces.  The  inclusion  of  electrodes  within  the  device  also  causes  a 
problem  due  to  the  non-linearity  of  the  equation  when  carrier  statistics  are 
included.  Finally,  the  method  solves  the  Poisson’s  equation,  but  cannot  be 
extended  to  the  Newton-step  form  of  the  equation  as  described  in  Chapter  3. 
The  linearization  term  resulting  from  the  Boltzmann  dependence  of  charge 
on  potential  re-introduces  the  spacial  coupling  of  the  equations  which  the 
FFT  method  eliminated.  It  is  possible  to  use  the  FFT  method  on  sub-regions 
of  the  device  and  then  couple  these  regions  together  using  an  iterative  tech¬ 
nique.  This  makes  it  possible  to  simulate  more  complex  device  structure's, 
but  some  of  the  advantages  of  a  direct  method  are  lost.  A  version  of  the  FFT 
method  is  used  in  the  CUPID  program  [1.4]  developed  by  J.  Frey. 

Of  the  several  iterative  techniques  for  solving  the  matrix  equations,  the 
most  commonly  used  is  successive  overrelaxation  (SOR)[4.5]  or  a  variant,  suc¬ 
cessive  block  overrelaxation  (S130R)  usually  implemented  as  successive  line 
overrelaxation  (SLOP).  These  methods  have  been  well  analyzed  in  the  litera¬ 
ture  [-1.6,  4.7],  so  they  will  be  discussed  only  briefly  here.  Theoretical  analysis 
of  the  SOR  and  SBOIl  methods  has  concentrated  on  rectangular  grids  and 
much  has  been  written  on  the  choice  of  optimum  relaxation  para  int  er  values. 
Although  the  SOR  techniques  are  applicable  to  any  type  of  grid,  the  optimum 
relaxation  parameter  value  is  not  easily  obtained  for  non-rectangulnr  grids, 
so  the  convergence  rate  of  the  method  suffers. 

Theoretically,  the  operation  count  of  SOR  is  0(/V3/2)  but  this  rate 
is  rarely  observed  in  practice  since  the  number  of  iterations  is  sensitive  to 
the  accuracy  of  the  initial  guess  of  the  solution  and  to  the  values  of  the 
coupling  coefficients.  The  convergence  rate  depends  on  the  eigenvalues  of 
an  iteration  matrix  derived  from  the  coefficient  matrix.  A  necessary  and 
sufficient  condition  for  convergence  is  that  the  largest  eigenvalue  must  be  less 
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than  one.  This  condition  is  guaranteed  if  the  coefficient  matrix  is  positive 
definite  (as  it  is  for  the  discretization  of  Poisson’s  equation  as  discussed  in 
Chapter  3.)  A  modified  SOR  method  is  used  in  the  TWIST  program  [4.8]. 

In  SOR,  the  solution  is  updated  for  each  grid  point  independently  and 
in  succession  based  on  current  values  of  the  solution  at  neighboring  points. 
In  SLOR,  the  solution  for  an  entire  line  of  points  is  computed  simultaneously 
based  on  the  values  of  the  solution  at  points  adjacent  to  this  line.  This 
has  the  advantage  of  providing  a  consistent  solution  to  the  equation  along 
a  line  of  points  instead  of  at  only  one  point.  TVs  a  result,  convergence  is 
generally  obtained  with  fewer  iterations.  The  disadvantage  is  that  one  must 
solve  a  set  of  simultaneous  equations  for  each  line.  The  disadvantage  is  not 
great,  however,  since  very  efficient  solution  techniques  exist  for  solving  the 
tridiagonal  matrix  equations  which  result.  Thus  there  is  generally  a  net 
decrease  in  the  amount  of  work  required  for  convergence  over  simple  SOR. 
Theoretically,  the  operation  count  is  O(Nrn)  where  m  is  the  number  of  lines 
in  the  iteration. 

Use  of  SLOR  normally  limits  one  to  a  rectangular  grid  since  the  iteration 
is  done  on  a  line  by  line  basis.  (There  is  nothing  to  prevent  selection  of  ser¬ 
pentine  “lines”  in  an  irregular  triangular  grid;  however,  there  is  a  penalty  in 
the  comput  ational  overhead  required  to  access  these  fi  les  and  the  convergence 
properties  of  such  a  method  are  unknown.)  For  IGFET  simulations  on  a  rec¬ 
tangular  grid,  the  lines  chosen  are  typically  vertical  lines  normal  to  the  gate 
since  the  potentials  are  more  strongly  varying  in  this  direction  than  in  the 
lateral  direction.  One  is,  in  effect,  solving  several  one-dimensional  problems 
side-by-side  and  coupling  them  together  with  the  iteration.  When  lateral 
potential  variations  become  large,  alternating  between  use  of  the  horizontal 
and  the  vertical  lines  may  provide  some  increase  in  convergence  rate;  however, 


the  increase  in  computational  overhead  offsets  most  of  the  advantage.  SLOR 
is  used  in  the  GEMINI  program  [4.9]. 

Another  iterative  method  commonly  used  in  semiconductor  device 
analysis  is  Stone’s  method  [4.10]  also  known  as  the  Strongly  Implicit 
Procedure  (SIP).  It  is  one  of  several  iterative  procedures  based  on 
modification  of  LU  decomposition  to  reduce  the  fill  problem  addressed  earlier. 
In  SIP,  an  L  and  U  matrix  are  computed  such  that  LU  =  M  +  E  where  both 
L  and  U  have  non-zero  elements  only  where  M  has  non-zero  elements,  thus 
there  is  no  fill.  The  elements  of  L  and  U  are  very  easily  computed  from  the 
elements  of  M.  The  matrix  E  has  the  same  non-zero  structure  as  A/  except 
for  two  additional  non-zero  diagonals.  The  combination  M  +  E  is  such  that 
the  iteration 

(M  +  E)xn+l  ={M  +  E)xn  -  (Mxn  -  y)  4.4 

converges  rapidly  to  the  proper  solution  where  i"  is  the  estimate  of  the 
solution  vector  at  the  nth  iteration. 

In  computing  the  elements  of  L  and  U ,  an  iteration  parameter  between 
zero  and  one  is  used  to  accelerate  the  convergence  of  the  algorithm.  The. 
choice  of  this  parameter  is  critical  since  too  large  a  value  will  result  in  diver¬ 
gence  and  too  small  a  value  will  result  in  slow  convergence.  Typically,  this 
paramter  is  cycled  through  a  range  of  values  in  order  to  insure  an  adequate 
degree  of  stability  while  maintaining  a  good  convergence  rate.  For  problems 
of  interest,  choosing  an  optimum  set  of  parameter  values  is  very  difficult 
so  conservative  estimates  are  made  yielding  sub-optimum  convergence.  The 
operation  count  for  this  method  is  0{N  In  N). 

The  method  was  developed  for  solving  equations  on  a  two-dimensional 
rectangular  grid;  however,  it  is  extendable  to  three-dimensional  grids  and  to 
multi-variable  problems.  The  principal  requirement  is  that  the  coefficient 
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matrix  have  a  highly  regular  banded  structure.  This  limits  the  grids  to 
rectangular,  rectangular  based  triangular,  or  cubic  structures.  Actually,  it 
is  the  highly  regular  connectivity  of  the  grid  that  is  important,  so  that  the 
grids  may  be  distorted  as  long  as  they  retain  their  regular  connectivity.  The 
SIP  method  is  used  in  the  CADDET  program  [4.11]. 

Another  method  based  on  modification  of  LU  decomposition  is  the 
Incomplete  Cholcsky — Conjugate  Gradient  (ICCG)  method  [4.12].  The  con¬ 
jugate  gradient  method  [4.13]  begins  with  an  estimate  of  the  solution  vec¬ 
tor  and  iteratively  improves  this  estimate  by  converging  to  the  exact  solu¬ 
tion  along  a  set  of  orthogonal  error  vectors.  Thus,  with  infinite  preci¬ 
sion  arithmetic,  the  exact  solution  would  be  obtained  in  N  iterations  since 
N  orthogonal  vectors  completely  span  the  solution  space.  In  fact,  if  the 
matrix  M  has  only  R  distinct  eigenvalues,  then  the  algorithm  must  con¬ 
verge  in  only  R  iterations.  Finite  precision  arithmetic  will  result  in  slower 
convergence,  but  since  exact  solutions  arc  generally  not  required,  sufficient 
accuracy  will  normally  be  achieved  in  the  number  of  iterations  described 
above.  Unfortunately,  the  coefficient  matrix  of  Poisson’s  equation  typically 
has  widely  spread  and  uniformly  distributed  eigenvalues  so  that  the  full  N 
iterations  are  required.  The  incomplete  Cholcsky  decomposition  algorithm, 
however,  provides  a  method  of  transforming  the  coefficient  matrix  into  an 
approximation  of  the  identity  matrix.  The  eigenvalues  of  this  approximate 
identity  matrix  are  highly  degenerate  so  only  a  few  are  distinct..  As  a  result, 
the  conjugate  gradient  method  converges  very  rapidly  when  used  with  the 
transformed  coefficient  matrix. 

For  a  symmetric,  positive  definite  matrix  (as  is  the  discretized  Poisson 
coefficient  matrix)  an  LLT  decomposition  may  be  obtained  instead  of  an  LU 
decomposition  where  L7  is  the  transpose  of  L.  This  is  the  Cholcsky  decom- 
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position,  LLt  =  M .  This  decomposition  results  in  the  same  fill  problem 
as  does  LU  decomposition;  however,  the  incomplete  Cholesky  decomposition 
avoids  this  fill  by  simply  zeroing  out  all  element  locations  in  L  which  are  zero 
in  M  so  that  L  retains  the  sparseness  of  M.  Since  the  elements  of  L  are 
computed  recursively,  the  zeroed  elements  influence  the  computation  of  later 
elements.  The  result  is  a  relation 

LLt  —  M  +  E  4.5 

much  like  that  obtained  for  SIP.  I’rc-  and  post-multiplying  by  inverses  yields 

L_1M(Lt)“1  =  I  -  L~lE{Lr)-1.  4.6 

If  the  error  matrix  term  on  the  right  hand  side  is  small,  then  the  left  hand 
side  is  approximately  the  identity  matrix  and  may  be  used  in  the  conjugate 
gradient  algorithm  for  rapid  convergence.  Physically,  M-1  (the  Greens 
function)  represents  the  coupling  between  a  node  and  its  neighbors.  Ignoring 
fill  in  the  L  matrix  is  effectively  neglecting  the  Greens  function  coupling  to 
the  distant  (non-adjacent)  neighbors  of  a  node. 

The  number  of  operations  per  iteration  of  this  method  is  comparable 
to  that  of  SBOR  or  SIP;  however,  the  rapid  convergence  (especially  for  stiff 
problems)  often  results  in  significantly  less  total  work.  The  method  does 
not  require  any  particular  matrix  structure  so  that  any  type  of  rectangular 
or  triangular  grid  may  be  used.  The  author  knows  of  no  device  simulation 
program  currently  employing  the  IGCG  method. 

4.1.2  Continuity  Equation 

The  continuity  equation,  as  discretized  in  Chapter  3,  does  not  yield  a 
symmetric  or  positive  definite  matrix.  As  a  result,  some  of  the  mat  rix  solution 
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algorithms  discussed  in  conjunction  with  the  solution  of  Poisson’s  equation 
cannot  be  applied  to  the  solution  of  the  continuity  equation.  Since  no  strong 
statement  may  be  made  concerning  the  eigenvalues  of  the  continuity  equation 
coefficient  matrix,  the  SOR  and  SBOR  methods  may  converge  very  slowly 
or  even  diverge  for  some  conditions.  The  FFT  method  is  also  not  applicable 
because  the  Scharfetter-Gummel  discretization  results  in  a  collection  of  node- 
to-node  differential  equations  for  the  transport  equation  rather  than  a  single 
global  differential  equation. 

LU  decomposition,  on  the  other  hand,  is  suitable  for  use  on  the  con¬ 
tinuity  equation  and  is  used  in  the  PISCES  program.  No  instabilities  have 
been  observed  in  using  this  method.  Both  the  SIP  and  ICCG  methods  also 
appear  suitable;  however,  a  modified  form  of  the  ICCG  method  is  required 
since  the  coefficient  matrix  is  not  symmetric.  For  asymmetric,  matrices,  the 
incomplete  Cholesky  decomposition  must  be  replaced  by  an  incomplete  LU 
decomposition,  where  the  incomplete  LU  decomposition  algorithm  is  fully 
analogous  to  that  for  the  incomplete  Cholesky  decomposition. 

4.1.3  Renumbering  Algorithms 

The  total  number  of  operations  and  data  storage  required  for  a  matrix 
solution  using  LU  decomposition  is  dependent  not  only  on  the  number  of 
nodes  in  the  grid  and  how  they  are  interconnected,  but  also  on  how  the 
nodes  are  numbered.  Consider,  for  example,  the  rectangular  and  rectangular- 
based  triangular  grids  of  Figures  2.01c  and  2.0 le.  If  the  nodes  in  these-  grids 
are  numbered  from  top  to  bottom  along  the  left  hand  column  of  w  nodes 
and  proceeding  column  by  column  from  left  to  right  along  the  ii  columns, 
then  the  structure  of  the  coefficient  matrices  will  appear  as  in  Figures  l.la 
and  4.11).  The  bandwidths  of  the  coefficient  matrices  ar(-  approximately  2m. 

(JO 
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Fig.  4.1.  Structure  of  coefficient  matrix  for:  (a)  rectangular  grid,  (b)  rectan¬ 
gular  based  triangular  grid.  The  clashed  lines  represent  partially 
filled  matrix  diagonals  resulting  from  the  choice  of  the  /  or  \  rec¬ 
tangle  diagonals  for  subdivision  into  triangles. 


Since  the  fill  associated  with  LU  decomposition  will  occur  totally  within  the 
bandwidth  of  the  matrix,  the  numbering  of  a  rectangularly  connected  grid 
should  proceed  in  the  direction  of  the  fewest  males,  by  rows  or  by  columns, 
in  order  to  minimize  the  bandwidth. 

Other  more  sophisticated  grid  numbering  schemes  exist  including  the 
method  of  nested  dissection  published  by  J.  A.  George  [4.14).  Using  this 
ordering,  the  operation  count  for  LU  decomposition  may  hi'  reduced  from 
0(N 2)  *o  0(/Vn)  and  the  needed  storage  from  O(A'n)  to  ()(A'ln/t).  An 
example  of  this  ordering  is  shown  for  a  7  by  7  rectangular  grid  in  Figure 
4.2.  Figure  4.2a  shows  the  partitioning  of  the  grid  all  nodes  labeled  1  are 
numbered  first  followed  by  those  labeled  2  and  then  3.  Figure  4.2b  shows  the 
a<  tual  numbering  with  the  resulting  matrix  structure  shown  in  Figure  4.3. 
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Fig.  4.2.  Nested  dissection  numbering  scheme  for  a  7  by  7  grid.  Part  (a) 
shows  the  partitioning  and  part  (b)  the  actual  numbering  of  (lie 
grid. 

Although  the  method  was  developed  for  rectangular  grids,  it  provides 
similar  benefits  for  the  rectangular-based  triangular  grids  used  in  PISCICS. 
Using  the  same  7  by  7  grid  and  numbering  as  in  Figure  4.2  but  subdividing 
the  grid  into  triangles,  the  matrix  structure  of  Figure  4.4  results.  The  struc¬ 
ture  is  essentially  the  same  as  that  of  Figure  4.3  and  results  in  comparable 
improvements  in  operation  count  and  storage. 

The  numbering  scheme  for  nested  dissection  works  perfectly  only  for 
square  grids  (i.e.  m  —  n )  and  only  when  m  —  2l  —  1  when'  l  is  an  integer, 
such  as  7  by  7,  15  by  15,  or  31  by  31  node  grids.  For  all  other  grids, 
slight  perturbations  in  the  numbering  are  required  in  order  to  ach'cve  the 
same  general  pattern  over  the  grid,  thus  there  are  many  possible  numbering 
patterns  for  any  given  grid.  Several  of  these  nested  dissection  numbering 
schemes  were  attempted  for  a  323  node  (17  by  ID)  PISCICS  grid.  Table  1.1 
shows  the  number  of  non-zero  entries  in  the  L  and  U  matrices  for  three  of 
these  numberings  plus  the  normal  column  numbering. 
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Fig.  4.3.  Coefficient  matrix  structure  for  the  nested  dissection  of  Figure  1.2 
on  a  rectangular  grid. 

The  symmetric,  numbering  is  a  perfectly  symmetric  dissection  about  the 
vertical  and  horizontal  centerlines  of  the  grid.  The  automatic  numbering  is 
roughly  the  same  but  not  as  symmetric.  It  is  generated  by  an  algorithm 
developed  by  I.  S.  Duff  [4.15]  and  modified  by  Prof.  R.  J.  Lomax  of  the 
University  of  Michigan.  The  truncated  numbering  is  simply  a  corner  of  a 
perfect  dissection;  that  is,  beginning  at  the  upper-left-hand  corner  node,  the 
perfect,  dissection  pattern  is  followed  horizontally  and  vertically  until  running 
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Fig.  4.4.  Coefficient  matrix  structure  for  the  nested  dissection  of  Figure  4.2 
on  a  rectangular  based  triangular  grid. 

out  of  grid.  The  results  show  that  the  symmetric  numbering  scheme  has 
the  greatest  reduction  in  storage  (20^)  followed  closely  by  the  automatic 
numbering.  The  truncated  numbering  results  in  an  increase  in  storage.  Note 
that  the  original  coefficient  matrix  contains  only  2223  non-zero  entries  while 
the  combined  L  and  U  matrices  of  the  best  numbering  contain  6732.  This 
difference  is  the  fill  generated  in  the  decomposition  process. 

A  decrease  in  storage  implies  a  decrease  in  the  number  of  operations 


Table  4.1 


Nested  Dissection  Efficiencies 


Numbering 

Scheme 

L+U  Storage 
(*) 

Solution  Time 
(sec) 

Normal 

8308 

2.0 

Symmetric  N.  D. 

6732 

1.8 

Automatic  N.  D. 

6948 

1.9 

Truncated  N.  D. 

8554 

2.9 

323  Node  grid  with  coefficient  matrix 
storage  =  2223  on  DEC-10  computer. 

required  since  fewer  fill  elements  arc  computed.  This  is,  in  fact,  the  case  as 
can  be  seen  in  the  seconds-pcr-solution  column  of  Table  1.1.  Note  that  this 
table  of  data  was  obtained  on  a  preliminary  version  of  PISCES  running  on 
a  large  time-share  computer  so  the  execution  times  arc  somewhat  biased  by 
load  variations.  Additional  timing  comparisons  were  obtained  on  a  525  node 
(21  by  25)  grid.  These  results  are  summarized  in  Table  1.2  and  show  that  the 
improvements  in  storage  and  solution  time  generally  agree;  with  that  predicted 
by  theory  and  r'so  that  the  improvements  arc  greater  for  larger  grids.  It 
appears  that  the  nested  dissection  method  provides  as  much  improvement 
for  the  rectangular  based  triangular  grids  of  PISCES  as  for  the  rectangular 
grids  for  which  it  was  developed. 

There  is  one  subtle  drawback  to  the  use  of  nested  dissection  for  moderate 
grid  sizes.  Table  4.3  shows  the  average  matrix  solution  time  for  PISCES  on 
the  I1P-1000E  for  three  different  conditions  scalar  arithmetic,  vector  arith¬ 
metic  and  vector  arithmetic  with  dissection.  The  PISCES  program  uses  a 


Table  4.2 


Nested  Dissection  Efficiencies  versus  Grid  Size 


L+U  Storage 

Solution  Time  ( 

sec) 

Grid  Size 

Normal 

Dissected 

Ratio 

I- 

Normal 

Dissected 

Ratio 

323 

8308 

6948 

1 

.84 

2.0 

1.9 

.96 

525 

20924 

12900 

.62 

5.8 

3.9 

.66 

Ratio 

2.5 

1.9 

2.9 

2.0 

Theory 

2.1 

1.8 

2.7 

2.1 

matrix  equation  solution  package  entitled  Vectorized  General  Sparsity  algo¬ 
rithms  (VEGES)  (4.16j  which  puts  all  vectorizable  steps  in  the  LU  decom¬ 
position  process  into  vector  operation  form.  On  true  vector  architecture  com¬ 
puters,  significant  time  savings  can  be  realized  by  exploiting  the  vectorizable 
operations;  however,  even  on  non- vector  machines  (such  as  the  HP-1000F) 
some  savings  can  be  achieved  with  microcoded  vector  instructions.  In  Table 
4.3,  the  scalar  version  me  program  uses  only  scalar  operations  while  the 
vector  version  uses  the  Vector  Instruction  Set  of  the  IUM000F  computer. 
The  vector  version  shows  a  factor  of  three  improvement  in  solution  time  over 
the  scalar  version.  When  nested  dissection  is  added  to  the  vectorized  version, 
however,  the  solution  time  increases.  This  increase  is  duo  to  a  reduction 
in  the  average  vector  length  caused  by  the  nested  dissection  renumbering. 
Depending  on  the  computer,  operations  on  short,  vectors  can  take  longer 
than  equivalent  scalar  operations  due  to  the  overhead  associated  with  vector 
operation  startup.  The  point  of  diminishing  returns  for  vector  versus  scalar 
operation  depends  solely  on  the  computer  being  used.  As  grid  size  increases, 
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Table  4.3 


Nested  Dissection  Effect  on  Vector  Arithmetic 


Conditions 

Solution  Time 
( sec ) 

Seal ar 

18 

Vector 

6 

Vector+ 

Dissection 

8 

HP-1000F  Computer 


however,  nested  dissection  provides  a  net,  solution  time  reduction  in  spite  of 
vector  shortening.  The  PISCES  program  contains  automatic  nested  dissec¬ 
tion  renumbering  as  a  user  controlled  option. 

Other  renumbering  schemes  were  also  attempted  including  various  forms 
of  diagonal  numbering  and  clustered  numbering.  None  of  these  methods 
showed  any  significant  reduction  in  storage;  in  fart,  most  showed  significant 
increases. 


4.2  Solution  of  Coupled  Equation s 


This  section  addresses  possible  ways  in  which  the  coupled  equations, 
Poisson  and  continuity,  can  be  solved  so  that  each  is  consistent  with  the  other. 
Factors  affecting  the  degree  of  coupling,  factors  influencing  the  convergence 
rates  of  iterative  methods  and  ways  of  accelerating  convergence  are  \>verod. 
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4.2.1  Solution  Methods 


Figure  4.5  shows  flow  charts  of  the  two  principal  approaches  to  solv¬ 
ing  the  coupled  equations —the  simultaneous  method  and  the  alternating  (or 
Gummels)  method.  The  comparison  between  these  two  methods  is  some¬ 
what  analogous  to  the  comparison  of  direct  versus  iterative  matrix  solution 
methods.  The  alternating  method  takes  less  work  per  pass  but  may  require 
many  passes  and  thus  more  work  for  a  consistent  solution  than  the  simul¬ 
taneous  method. 

The  simultaneous  method  involves  appending  the  two  matrix  equations 
together  to  form  a  single  matrix  equation  which  is  twice  as  large.1  In  addition, 
the  partial  derivatives  of  all  combinations  of  potential  and  carrier  concentra¬ 
tion  for  adjacent  nodes  must  he  included.  This  doubles  the  number  of  rows 
and  approximately  doubles  the  number  of  non-zero  entries  per  row  of  the 
matrix  since  the  carrier  concentration  at  each  node  is  coupled  not  only  to 
the  carrier  concentrations  of  each  adjacent  node,  but  also  to  the  potentials. 
The  result  is  a  matrix  equation  with  nearly  four  times  the  number  of  non¬ 
zero  entries  as  either  coefficient  matrix  alone,  and  more  than  four  times  as 
many  operations  per  solution.  This  rapid  multiplication  of  elfort  is  prohibi¬ 
tive  for  grids  of  any  practical  size.  Additionally,  the  compulation  of  the  par¬ 
tial  derivatives  becomes  more  difficult  as  higher  order  physical  phenomena 
such  as  field  dependent  mobility  are  introduced.  For  these  reasons,  most 
device  simulation  program  designers  have  shunned  the  simultaneous  method; 
however,  some  insist  that  simultaneous  solutions  are  imperative  for  certain 
bias  conditions  [4.17]. 

1  For  two-carrier  simulation,  both  continuity  equation  matrices  would  be  ap¬ 
pended  resulting  in  a  matrix  three  times  as  large. 
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Fig.  •1.5.  Algorithm  (low  of  simultaneous  and  alternating  methods  for  solu¬ 
tion  of  the  coupled  equations. 

The  alternating  method  of  Gummel  [-1 .  IS]  is  the  most  commonly  used 
method  for  obtaining  consistent  solutions  to  the  Poisson  and  continuity  equa¬ 
tions.  As  illustrated  in  the  flow  chart  of  Figure  1.5,  beginning  with  an  initial 
guess  of  potentials  and  carrier  concentrat  ions,  the  Poisson  equation  is  sohed 
for  node  potentials.  The  carrier  concentrations  arc  updated  based  on  the  new 
potentials  and  appropriate  carrier  statistics  (Holt/mann  or  Fermi-Dirac)  and 
the  Poisson  equation  is  solved  again.  This  procedure  is  repeated  until  the 
potential  change  is  below  some  convergence  criteria  limit.  This  is  the  inner 

78 


loop  of  the  How  chart  of  Figure  4.5b. 

Implicit  in  the  updating  of  carrier  concentrations  based  on  potential 
changes  and  carrier  statistics  is  the  assumption  of  a  fixed  quasi-Fermi  level. 
In  fact,  the  quasi-Fermi  level  is  an  unknown  which  must  be  determined. 
This  is  done  implicitly  in  the  outer  loop  by  solving  the  continuity  equation 
and  updating  the  carrier  concentrations  without  changing  the  potentials. 
These  new  carrier  concentrations  require  that  the  Poisson  equation  be  solved 
again  so  the  algorithm  loops  back  to  the  top  of  the  flow  chart.  Eventually, 
potentials  and  carrier  concentrations  converge  and  are  consistent  with  both 
the  Poisson  and  continuity  equations. 

The  principal  advantage  of  the  alternating  method  is  that  less  work 
is  expended  on  each  pass  through  the  outer  loop  than  in  the  simultaneous 
method.  A  significant  disadvantage,  however,  is  that  the  convergence  rate 
of  the  method  is  dependent  on  the  device  operating  conditions  and  may  be. 
very  slow.  In  the  simulations  run  on  the  PISCES  program,  for  example,  the 
number  of  outer  loop  iterations  required  for  convergence  without  acceleration 
varied  from  one  to  sixty  depending  on  the  device  biasing  conditions.  The 
reduction  of  this  large  number  of  iterations  was  a  major  thrust  of  the  present 
work. 

4.2.2  Convergence  Acceleration  for  the  Alternating  Method 

Table  4.4  shows  a  typical  convergence  pattern  for  a  MOSFET  biased 
below  threshold.  The  left,  hand  column  is  the  iteration  count  of  the  outer 
loop  of  the  alternating  method,  the  center  column  is  the  error  measure  for 
each  iteration  of  the  inner  loop  (Newton  iteration  on  Poisson’s  equat  ion),  and 
the  right  hand  column  is  the  error  measure  for  each  iteration  of  the  outer 
loop.  The  inner  loop  error  measure  is  the  average  of  the  absolute  values  (the 


79 


Table  4.4 


Subthreshold  Convergence 


Outer 

Loop 

Count 

Inner 

Loop 

Error 

_ 

Outer 

Loop 

Error 

3 

7  E  -  4 

5 

8E-12 

1 

3 

7E- 

4 

1 

5  E  - 1 

1 

4E-2 

1 

7E-3 

3 

4E-4 

1 

3E-5 

3 

3E-8 

2 

1 

4E- 

2 

1 

2E-12 

3 

1 

2E- 

12 

Total 

9 

3 

Total  matrix  solutions  =  12 
Gate  =  .5V,  Drain  =  ,01V 


one  norm  of  numerical  analysis)  of  the  incremental  potentials  resulting  from 
the  Poisson  solution.  The  outer  loop  error  measure  is  the  one  norm  of  the 
net  change  in  the  node  potentials  from  one  pass  through  the  outer  loop  to 
the  next. 

Analyzing  the  convergence  of  Table  4.4,  on  the  first  pass  through  the 
outer  loop,  Poisson’s  equation  converges  in  two  iterations.  The  continuity 
equation  solution,  however,  alters  the  carrier  concentrations  and  results  in 
six  Poisson  iterations  on  the  second  pass  through  the  outer  loop.  On  the 
third  pass,  only  one  Poisson  solution  is  required  and  the  change  in  potential  is 
extremely  small  indicating  that  the  algorithm  has  converged.  Twelve  matrix 
solutions  are  required  for  this  simulation,  nine  for  Poisson  and  three  for 
continuity.  This  rapid  convergence  is  due  to  the  very  weak  coupling  between 


80 


the  Poisson  and  continuity  equations.  When  a  MOSFET  is  biased  below 
threshold,  the  dominant  charge  (besides  boundary  charge)  is  the  space  charge 
of  ionized  impurities  in  the  depletion  regions.  Since  this  charge  is  immobile, 
the  continuity  equation  has  little  effect  and  a  Poisson  solution  is  essentially 
all  that  is  required. 

Another  point  worth  noting  in  Table  4.4  is  the  convergence  rate  of  the 
inner  loop.  Newton’s  method  has  quadratic  convergence  which  means  that 
when  the  solution  estimate  is  in  the  vicinity  of  the  correct  solution,  the  error 
in  the  solution  estimate  will  be  squared  with  each  iteration.  This  type  of 
behavior  is  evident  in  the  inner  loop  error  measure  shown  in  the  table. 

A  different  convergence  pattern  is  seen  in  the  MOSFET  linear  region 
simulation  of  Table  4.5.  In  this  case,  there  is  an  inversion  layer  of  free 
carriers  at  the  semiconductor  surface  resulting  in  stiff  coupling  between  the 
Poisson  and  continuity  equations.  Only  pieces  of  the  complete  convergence 
sequence  arc  shown,  but  it  is  clear  that  each  continuity  solution  alters  the 
carrier  concentrations  enough  to  negate  the  previous  Poisson  solution  and 
cause  several  more  iterations  of  the  inner  loop.  A  total  of  151  inner  loop  and 
46  outer  loop  iterations  are  required.  It  is  this  behavior  which  undermines 
the  utility  of  the  alternating  method  in  device  simulation.  In  the  search  for 
ways  to  reduce  the  simulation  time  required  for  devices  above  threshold,  four 
technique,  were  derived  which  may  be  used  to  reduce  this  time  by  a  factor 
of  three  to  four.  These  techniques  are  detailed  in  the  following  paragraphs. 

4.2.2. 1  Projection  of  the  Initial  Guess 

The  first  convergence  acceleration  technique  is  the  use  of  projection  or 
the  initial  guess  from  previous  solutions.  In  t he-  PISCES  program,  the  first 
bias  condition  solved  for  newly  generated  device  grids  is  the  Hat  band  case. 
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Table  4.5 


Linear  Region  Convergence 


Outer 

Inner 

Outer 

Loop 

Loop 

Loop 

Count 

Error 

Error 

1.0E-2 

2.GE-3 

2.  IE -3 
6.1E-4 

1. IE-4 
4.2E-6 

1 

4.5E-2 

2.6E-2 

1.3E-2 

2.3E-3 

1. IE-3 
3.5E-4 
4.3E-5 

6.9E-3 

2 

. 

2.2E-2 

16 

* 

* 

3.3E-3 

6.6E-4 

6. 6  E - 5 

2.8E-3 

17 

2.5E-3 

• 

23 

1.0E-3 

9.0E-5 

1 . IE-3 

24 

* 

9. 4E-4 

45 

1.0E-4 

8.3E-5 

46 

8.3E-5 

Total 

151 

— _ — _ 

46 

Total  matrix  solutions  =  197 
Gate  =  2.7V,  Drain  =  1 . OV 
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The  initial  guess  for  this  case  is  the  charge  neutral  condition.  When  not  using 
the  projection  technique,  the  initial  guess  for  every  simulation  thereafter  is 
the  previous  solution.  If  a  sequence  of  bias  steps  is  being  simulated,  for 
example,  each  solution  is  used  as  the  initial  guess  for  the  next  bias.  The 
result  is  slow  initial  convergence  and  occasional  instability  if  too  large  a  bias 
step  is  attempted. 

The  projection  method  involves  extrapolating  the  initial  guess  for  the 
new  bias  condition  based  on  the  solutions  at  two  previous  bias  conditions. 
Only  one  contact  bias  may  be  varying  between  the  two  previous  biases  and 
the  new  bias.  The  assumption  is  that  the  potentials  and  quasi-F<  rmi  levels 
for  each  node  will  vary  linearly  with  the  bias.  Given  two  solution  files  for 
biases  of  Vx  and  V2  and  a  new  bias  of  V3,  an  extrapolation  factor  is  defined 
as 


a 


V3  -  V2 
V2  Vi 


(-1.7) 


The  projected  initial  guess  for  potential  and  quasi-Fermi  level  at  the  new  bias, 
V’3  and  03,  is  then  computed  for  every  node  in  the  device  by  the  relations 


03  —  02  +  «(02  -  01 )  (TK) 

and 

03  ---  02  +  "(02  01  )•  ('!•■!) 


It  is  easil)  shown  that  these  relations  lead  to  an  extrapolation  of  the  electron 
concentration  (assuming  Holtzmann  statistics)  of 


n3  —  «a 


(-1.10) 


These  relations  provide  very  good  initial  guesses.  In  charge  neutral 
regions,  the  estimates  are  very  accurate.  In  depletion  regions  the  estimates 
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arc  not  as  good,  but  since  there  is  little  mobile  charge  there  the  coupling 
between  Poisson’s  equation  and  the  continuity  equation  is  weak  so  that  those 
errors  which  do  exist  are  quickly  corrected  by  Poisson  solutions.  Only  at 
depletion  edges  and  inversion  layers  do  the  estimates  produce  significant 
errors,  but  even  there  the  estimates  provide  a  smooth  stable  initial  guess 
with  excellent  convergence  probability.  The  combination  of  these  properties 
serve  not  only  to  speed  convergence  but  also  to  greatly  increase  the  allowable 
bias  step  size.  Another  advantage  of  this  method  is  its  simplicity  since  all 
regions  of  a  device  are  processed  identically. 

The  projection  method  provides  the  greatest  convergence  acceleration 
when  stepping  the  drain  bias  of  devices  in  saturation,  typically  reducing  the 
number  of  iterations  by  a  factor  of  two.  It  is  nearly  as  effective  in  the  linear 
region  but  tends  to  lose  some  of  its  effectiveness  in  projecting  initial  guesses 
through  the  transition  from  linear  to  saturation.  For  subthreshold  bias 
conditions  the  method  does  not  significantly  reduce  the  number  of  Herat  ions; 
however  it  does  allow  larger  bias  steps  to  be  taken  without  loss  of  stable 
convergence. 

4. 2. 2. 2  Single  Poisson’s  Equation  Iteration 

The  second  convergence  acceleration  technique  is  the  use  of  a  single 
Poisson  solution  per  outer  loop  iteration.  An  inspection  of  Table  1.5  reveals 
that  each  continuity  equation  solution  severely  alters  the  previous  Poisson 
solution  as  indicated  by  the  large  inner  loop  error  seen  after  each  outer  loop 
iteration.  In  short,  the  accurate  Poisson  solution  achieved  through  several 
iterations  of  the  inner  loop  is  unnecessary.  Hy  performing  only  one  inner 
loop  iteration  on  each  pass  through  the  outer  loop,  the  number  of  outer  loop 
iterations  increases  by  roughly  20f7  but  the  total  number  of  matrix  solutions 
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Fig.  4.6.  Algorithm  flow  of  the  single  Poisson  acceleration  scheme  showing 
elimination  of  the  inner  loop. 

decreases  by  roughly  40%.  Figure  1.6  shows  the  revised  flow  chart  with  the 
elimination  of  the  inner  loop.  This  technique  applies  only  to  the  linear  and 
saturation  bias  conditions.  Its  use  in  suhthreshold  simulations  will  generally 
increase  the  total  number  of  matrix  solutions  required. 

The  success  and  stability  of  this  met  hod  provokes  the  quest  ion  of  whether 
this  merging  of  the  two  iterative  loops  could  be  taken  a  step  further  when 
iterative  methods  are  used  for  the  matrix  equation  solutions.  That  is,  rather 
than  fully  converging  on  an  iterative  solution  to  the  Poisson  matrix  equation 
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and  then  fully  converging  on  a  solution  to  the  continuity  matrix  equation, 
perhaps  the  program  should  alternate  between  a  few  iterations  of  each  matrix 
solution.  This  possibility  is  worthy  of  additional  consideration,  but  was  not 
pursued  in  this  work. 

4. 2. 2. 3  Overrelaxation 

The  third  convergence  acceleration  technique  is  the  use  of  a  form  of  over- 
relaxation.  This  technique  must  be  used  only  with  the  single  Poisson  iteration 
described  above.  The  method  was  developed  by  observing  the  details  of  the 
convergence  of  simulations  using  projection  and  a  single  Poisson  solution. 
Figure  1.7  shows  the  potential  and  quasi-Fermi  potential  convergence  of  a 
node  in  the  channel  region  of  a  MOSFFT  biased  in  saturation.  The  error  is 
plotted  versus  the  outer  loop  iteration  count,  thus  each  iteration  represents 
two  matrix  solutions.  Finer  detail  of  the  potential  and  quasi-Fermi  potential 
convergence  is  shown  in  Figure  4.8  along  with  the  electron  concentration  for 
the  same  node  and  the  total  drain  current.  Note  that  the  latter  two  lag  by 
about  20  iterations. 

The  monotonic  nature  of  the  potential  convergence  for  this  node  over  t lie 
first  10  to  20  iterations  is  typical  of  nodes  observed  in  other  regions  of  the 
device.  Since  the  potential  increments  are  nearly  constant  for  every  iteration, 
it  appears  that  faster  convergence  can  be  obtained  by  merely  increasing  the 
size  of  tin'  potential  increments.  This  is  somewhat  analogous  to  overrelaxation 
in  iterative  matrix  solution  methods.  When  the  vector  of  potential  increments 
out  of  the  Poisson  solution  is  multiplied  by  a  factor  greater  than  one  before 
being  added  to  the  previous  potentials,  faster  convergence  does  indeed  result. 

Figure  4.9  shows  the  improvement  in  drain  current  convergence  obtained 
using  overrelaxat ion.  A  factor  of  1.0  means  no  overrelaxation  and  1.5  or  1.9 
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Fig.  4.7.  Potential  and  quasi-Fermi  potential  convergence  of  a  node  in  tin 
channel  region  of  a  MOSFET  biased  in  saturation. 
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Fig.  4.8.  Finer  detail  of  the  convergence  shown  in  Figure  4.7  including  node 
electron  concentration  convergence  and  total  drain  current  corner- 
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Fig.  4.9. 


Drain  current  convergence  acceleration  using  overrelaxation. 


mean  that  the  Poisson  solution  vector  is  multiplied  by  1.5  or  1.9  each  time 
before  updating  the  potentials.  Factors  of  two  or  larger  cause  immediate 
instability  in  the  iterations,  suggesting  a  similarity  between  this  method  and 
the  successive  overrelaxation  matrix  iterative  technique. 

The  PISCES  program  uses  a  variable  overrelaxation  factor  which  changes 
every  three  iterations,  typically  starting  at  a  large  value  and  then  decreas¬ 
ing.  This  factor  is  computed  using  an  algorithm  published  by  Carre  [1.19]; 
however,  fixed  factors  between  1.5  and  1.9  generally  work  equally  as  well. 
The  overrelaxation  technique  can  be  used  only  when  using  the  single  Poisson 
method  also,  thus  it  is  applicable  only  to  the  linear  and  saturation  regions  of 
operation. 


4.2.2.-1  Linearization  Term  Reduction 

The  final  convergence  acceleration  technique  is  reduction  of  the  lineariza¬ 
tion  term  in  Poisson’s  equation.  This  method  is  closely  related  to  the  over- 
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relaxation  method  and  their  use  is  mutually  exclusive.  The  basis  for  this 
method  can  be  derived  by  inspecting  Figure  4.7  and  observing  that  the  quasi- 
Fermi  potential  closely  tracks  the  potential  as  they  both  converge.  After  each 
update  of  the  potential  by  a  Poisson  solution,  the  continuity  equation  results 
in  a  nearly  identical  step  in  the  quasi-Fermi  potential.  This  violates  one  of 
the  assumptions  made  in  the  discretization  of  Poisson’s  equation  performed 
in  Chapter  3.  In  that  discretization,  the  linearization  of  the  carrier  statistics 
was  based  on  the  assumption  that  the  quasi-Fermi  level  would  remain  con¬ 
stant  as  the  potential  changed  resulting  in  a  predictable  change  in  electron 
concentration.  W  hen  using  the  single  Poisson  method,  however,  the  quasi- 
Fermi  potential  follows  the  potential  change  at  each  iteration;  as  a  result  the 
linearizing  term  in  the  discretization  should  he  made  smaller. 

The  Boltzmann  carrier  statistics  assumed  for  this  work  are  expressed  for 
electrons  as 

n  =  n,e(v'_‘>n)/v'T.  (4.11) 

The  linearization  of  Poisson’s  equation  then  uses  the  partial  derivative 


dn  n 
dd  =  VT  ’ 


(1.12) 


however,  if  Qn  tracks  d,  then  in  one  iteration  of  the  outer  loop, 


An  n 
A  ip  ^  Vt 


(4.13) 


Thus,  more  rapid  convergence  should  be  obtained  if  one  uses  a  linearizat ioi 


dn  n 

dd  =  °  VV’ 


(4.14) 


where  the  accelerating  factor,  a,  has  been  empirically  determined  to  lie  in 
the  range  .3  <  a  <  .6. 
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This  technique  is  roughly  equivalent  to  the  overrel notation  technique 
except  that  its  effects  are  limited  to  regions  of  the  device  with  high  electron 
concentration.  In  low  electron  concentration  regions,  the  linearization  term 
is  small  and  the  Poisson  solution  values  are  unperturbed.  In  high  electron 
concentration  regions.  the  linearization  term  dominates  and  the  size  of  the 
update  potential  solution  is  increased  by  a  factor  of  approximately  l/o. 

Figure  1.10  shows  the  drain  current  convergence  acceleration  achieved 
with  this  method.  The  normal  unacceleratcd  case  is  obtained  with  a  factor  of 
1.0.  A  factor  of  .5  is  comparable  to  selective  application  of  an  overrelaxation 
factor  of  2.0  and  is  seen  to  reduce  the  number  of  iterations  by  about  one 
half.  The  third  example  shown  used  a  factor  which  started  at  .2  for  the  first 
iteration,  and  increased  by  .0.33  each  iteration  until  reaching  .6  where  it  was 
fixed.  The  initial  convergence  was  very  rapid  but  unstable  causing  a  residual 
ringing  in  the  drain  current  even  after  the  acceleration  factor  had  moved 
into  the  stable  range  (>.5)  on  the  tenth  iteration.  The  convergence  rate  in 
this  case  was  only  marginally  better  than  that  obtained  using  a  fixed  value 
of  .5,  but  the  rapid  initial  convergence  obtainable  from  use  of  an  unstable 
acceleration  factor  value  will  generally  result  in  a  savings  of  a  few  iterations. 
An  acceleration  factor  starting  at  .3  and  increasing  by. 01  each  iteration  until 
reaching  .0  is  recommended. 

As  mentioned  earlier,  this  method  may  bo  used  interchangeably  with 
the  overrelaxation  method,  but.  they  cannot  be  used  simultaneously.  The 
performance  of  the  two  methods  is  nearly  identical  and  no  recommendations 
are  made  as  to  the  use  of  one  over  the  other. 

It  is  important  to  note  that  Figures  <1.9  and  1.10  show  the  drain  current 
error  within  a  range  of  plus  or  minus  ,\%.  Drain  current  accuracies  required 
for  device  simulation  are  more  typically  1 .0*  (  and  certainly  not  less  then 
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Fig. '1.10.  Drain  current  convergence  acceleration  using  reduction  of  the 
linearizing  term. 

0.1%.  Looking  at  the  figures,  these  levels  oT  convergence  are  achieved  in  12 
and  25  iterations  respectively,  thus  the  convergence  to  useful  levels  is  quite 
rapid  when  using  the  convergence  acceleration  techniques. 

The  convergence  acceleration  obtainable  from  combinations  of  these  ac¬ 
celeration  methods  on  a  MOSF1CT  operating  in  the  saturation  mode  is  ex¬ 
hibited  in  Table  4.6.  The  left  hand  column  shows  the  five  drain  biases  at 
which  the  simulations  were  run.  The  second  column  shows  the  number  of 
matrix  solutions  (Poisson  plus  continuity)  required  for  convergence  at  each 
drain  bias  with  no  acceleration.  The  total  number  required  for  all  five  biases 
is  shown  at  the  bottom  along  with  a  ratio  of  the  total  number  f’oi  each  column 
to  the  total  of  the  unaccelerated  column.  The  third  column  shows  the  ac¬ 
celeration  achieved  using  projection  alone.  The  fourth  and  fifth  columns  show 
the  acceleration  achieved  using  the  single  Poisson  iterations  alone  and  with 
projection.  The  final  column  shows  the  acceleration  achieved  by  using  all 
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Table  4.6 


Number  of  Matrix  Solutions  to  Convergence 


Acceleration 

Method 

Drain 

Bias 

None 

Projection 

1-Poisson 

1-Poisson 
Project i on 

Overrel ax  . 

1-Poisson 
Project  ion 

1.5 

189 

98 

120 

74 

68 

2.0 

188 

110 

86 

84 

68 

2.5 

183 

103 

88 

84 

48 

3.0 

183 

78 

88 

70 

58 

3.5 

184 

82 

90 

52 

48 

Total 

927 

471 

472 

364 

290 

Ratio 

1.00 

.51 

.51 

.39 

.31 

three:  projection,  single  Poisson,  anti  overrelaxation.  Comparable  data  was 
not  taken  for  the  linearizing  term  reduction  method,  but  results  are  quite 
similar  to  the  overrelaxation  results. 

The  data  in  this  table  is  slightly  distorted  in  the  conservative  direction 
due  to  the  convergence  criteria  which  was  being  used  at  the  time  it  was 
gathered.  The  iterations  were  stopped  when  the  change  in  drain  current 
from  one  iteration  to  the  next  became  relatively  small.  This  was  determined 
later  to  be  a  poor  criteria  since  it  usually  stopped  on  the  peak  of  the  drain 
current  overshoot  a  local  maximum  of  drain  current  error  magnitude  as  seen 
in  Figures  4.6  and  4.10.  Occasionally,  however,  the  algorithm  missed  the 
peak  and  stopped  many  iterations  later  resulting  in  a  more  tightly  converged 
solution.  Typically,  the  algorithm  stopped  on  the  peak  for  slower  converging 
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methods  but  not  on  faster  ones,  thus  the  entries  in  the  unacceleratcd  column 
of  Table  4.6  are  loosely  converged  while  some  of  the  entries  in  the  remaining 
columns  are  tightly  converged.  As  a  result,  the  net  improvement  ratio  is 
underestimated.  The  outer  loop  error  norm  as  described  earlier  is  a  more 
stable  measure  of  convergence  and  is  used  in  the  present  version  of  PISCES. 
A  convergence  criteria  for  this  error  norm  of  10-4  results  in  loose  convergence 
and  10-5  in  tight  convergence.  Use  of  this  convergence  criteria  for  all  of 
the  data  in  Table  4.6  would  show  a  total  acceleration  improvement  ratio  of 
approximately  .25. 

4.2.3  Convergence  Rate  Sensitivity  to  Bias  Conditions 

The  previous  discussions  have  referred  to  the  sensitivity  of  convergence 
rate  to  the  device  region  of  operation.  Figures  1.11  and  4.12  show  this 
sensitivity  in  the  computation  of  lp/Vc  and  lp/Vps  characteristics  for  a 
short  channel  ICFET  using  a  700  node  grid.  In  both  figures  the  drain 
current  is  shown  as  a  solid  line  and  the  total  solution  time  per  bias  point  as 
dots.  The  solution  time  was  measured  on  an  HP- 10001'’  minicomputer  using 
vector  arithmetic  with  each  matrix  solution  of  the  700  node  grid  requiring 
approximately  30  seconds.  Note  that  1  he  piecewise  nat  ure  of  1  he  drain  current 
curves  results  from  the  large  discrete  bias  steps. 

Figure  4.11  shows  the  drain  current  for  a  drain  bias  af  .01 V  and  a  gate 
bias  at  steps  of  ,2V  from  sub-threshold  to  the  linear  region  of  operation. 
Projection  of  the  initial  guess  was  the  only  convergence  acceleration  procedure 
used  for  these  simulations.  The  very  small  drain  bias  in  these  simulations 
avoids  the  need  for  any  other  acceleration  procedures.  The  solution  times 
are  small  for  all  bias  points  but  increase  slightly  near  threshold  where  the 
device  is  changing  from  subthreshold  to  linear  operating  characteristics.  This 
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increase  is  due  to  the  errors  incurred  in  projecting  an  initial  guess  in  a  region 
where  the  device  is  changing  modes  of  operation. 

Figure  4.12  shows  data  for  the  same  device  with  a  gate  bias  of  2V 
(well  above  threshold)  and  drain  bias  at  .5V  steps.  The  projection,  single 
Poisson,  and  linearizing  term  reduction  methods  of  convergence  acceleration 
were  used  for  these  simulations.  The  noticeable  aspects  of  this  data  are  the 
large  solution  times  required  for  drain  biases  at  the  knee  of  the  curve  and 
the  decrease  in  solution  time  as  the  drain  bias  increases.  The  first  data 
point  shown  on  this  figure  is  ,5V  and  occurs  in  the  transition  region  between 
the  linear  and  saturation  modes  of  operation  of  the  device,  thus  the  lower 
solution  times  required  in  the  linear  region  are  not  shown.  A  small  part  of 
the  reduction  in  solution  times  as  drain  bias  increases  is  the  improvement  in 
accuracy  of  the  projected  initial  guess,  but  the  iteration  output  reveals  that 
this  accounts  for  only  a  few  iterations  difference  per  bias  point.  The  principal 
difference  is  the  rate  of  convergence.  In  the  transition  region  it  takes  more 
than  ten  iterations  to  reduce  the  error  norm  by  one  order  of  magnitude,  but 
in  the  saturation  region  it  takes  less  than  five.  The  cause  of  this  performance 
is  unknown;  although,  it  appears  to  be  a  characteristic  of  the  convergence 
acceleration  schemes  since  it  is  not  observed  in  unaccelerated  simulations. 

Some  added  insight  is  provided  by  Figure  4.K1  which  shows  the  error  in 
surface  potential  at  each  of  the  first  IS  iterations  for  a  device  with  a  1.5  /nil 
gate  length.  Only  the  projection  and  single  Poisson  acceleration  schemes 
were  used  in  this  simulation.  It  is  evident  that  the  error  is  dominated  by  the 
first  harmonic  in  spaeial  frequency  between  the  source  and  drain  and  that  it 
converges  very  slowly.  Apparently  the  alternating  method  provides  adequate 
local  coupling  (adjacent  nodes)  for  the  Poisson  and  continuity  equations  but 
is  less  adequate  globally,  specifically  from  the  source  to  the  drain  along 
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the  inversion  channel.  It  would  appear  prudent  to  investigate  methods  of 
supplementing  this  coupling  across  the  length  of  the  channel. 

The  approximate  solution  times  for  the  various  device  operating  regions 
are  summarized  in  Table  '1.7.  An  average  number  of  matrix  solutions  required 
for  convergence  in  each  region  is  given  along  with  the  number  of  minutes 
required  for  three  different  grid  sizes. 

«H> 


Table  4.7 


Approximate  Solution  Times  on  IIP-1000F 


Operating 

Region 

Number 
of  Matrix 
Solutions 

Subthreshold 

7 

Linear 

25 

Transition  | 

50 

Deep  Saturation 

■ 

15 

Minutes  per  Bias  Point 


9.6  12.5 

19.1  25.0 


4.3  Summary 


A  variety  of  direct  and  iterative  matrix  solution  methods  are  presented 
and  their  applicability  to  the  Poisson  and  continuity  equations  arc  discussed. 
The  advantages  obtained  from  various  grid  renumbering  schemes  are  also 
discussed.  The  nested  dissection  renumbering  is  shown  to  provide  significant 
savings  in  both  storage  and  solution  time,  especially  for  largei  grids.  These 
improvements  are  tempered  by  the  fact  that  the  re-ordering  results  in  shorter 
vectors  in  the  solution  algorithm.  The  shorter  vectors,  in  turn,  slow  down 
the  processing  speeds  available  on  modern  vector  computers. 

Solution  of  the  coupled  system  of  equations  is  addressed  and  the  tradeoffs 
of  simultaneous  versus  alternating  solutions  are  discussed.  Convergence  of 
the  alternating  method  without  acceleration  is  analyzed  and  shown  to  hr 
prohibitively  slow  Tor  lGFPTs  biased  in  the  linear  and  saturation  modes.  Four 
convergence  acceleration  techniques  are  presented  which,  in  combination, 
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increase  the  convergence  rate  by  roughly  a  factor  of  four.  These  methods 
involve  computation  of  an  improved  initial  guess,  elimination  of  excessive 
solutions  of  Poisson’s  equation,  overrelaxation  of  the  potential  updates,  and 
reduction  of  the  Poisson  linearization  term. 

The  convergence  rate  is  examined  as  a  function  of  bias  for  operation 
below  threshold  and  operation  above  threshold  with  acceleration.  Solutions 
above  threshold  take  roughly  two  to  six  times  as  long  as  subthreshold  solu¬ 
tions.  Convergence  is  the  slowest  for  devices  biased  in  the  transition  region 
between  linear  and  saturation  operation  but  improves  as  the  device  is  biased 
deeper  into  saturation.  The  dominant  electrostatic  potential  error  observed 
during  convergence  above  threshold  is  a  first  harmonic  in  spacial  frequency 
between  the  source  and  drain  at  the  surface. 

The  next  chapter  provides  two  application  examples  of  the  PISCPS 
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Chapter  5 


APPLICATIONS 

This  chapter  covers  two  examples  of  applications  of  the  PISCES  pro¬ 
gram.  In  the  first  example,  an  NMOS  transistor  with  a  channel  implant 
is  simulated  in  the  punchthrough  region  of  operation  using  three  different 
simulation  programs —PISCES,  GEMINI  [5- 1  ] ,  and  CADDET  jo. 2].  The 
effects  on  punchthrough  current  of  varying  the  source/drain  junction  depths 
arc  also  shown.  The  second  example  demonstrates  the  utility  of  the  PISCES 
program  in  evaluating  physical  models.  The  field  dependent  mobility  of  field 
effect  devices  is  explored  through  use  of  a  distance-from-the-surface  mobility 
variation. 


5.1  MOSFET  Punchthrough 

A  potential  contour  plot  of  the  N-channel  MOSFET  simulated  for  this 
comparison  is  shown  in  Figure  5.1.  The  important  device  parameters  are: 


simulation  length 

3.5  pm 

simulation  depth 

5  pm 

oxide  thickness 

500  A 

channel  length 

.7  pm 

substrate  doping 

2  X  1015 

junction  depth 

A  pm 

channel  implant  dose 

3,1  X  10 

channel  implant  depth 

.15  pm. 
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Fig.  5.1.  Equi potential  contour  plot  of  N-channel  MOSFET  with  implanted 
channel. 

The  planar  oxide  was  used  in  order  to  accomodate  the  CADDET  program 
comparison. 

Figure  5.2  shows  the  results  of  the  punchthrough  simulations  using 
PISCES,  GEMINI,  and  CADDET.  The  drain  bias  was  varied  in  one  volt  steps 
from  one  to  ten  volts  while  holding  the  gate  at  —  .5VE  An  analysis  of  equi- 
potential  contour  plots  would  reveal  that  the  punchthrough  current  path  is  at 
the  surface  for  Vps  <  7E  and  in  the  bulk  for  Vps  >  7V.  The  A  fim  source 
and  drain  regions  are  modeled  as  Gaussian  implants  with  the  concentration 
peak  exactly  at  the  semiconductor  surface.  The  .45/im  junction  dept  h  device 
source  and  drain  have  the  same  Gaussian  shape  but  the  peak  is  shifted  .05 /un 
below  the  semiconductor  surface.  Only  PISCES  was  used  to  simulate  the 
.45  fim  junction  device.  The  punchthrough  current  for  this  device  exceeds 
that  of  the  .4  /un  device  by  as  much  as  a  factor  of  50.  This  fifty-fold  in- 
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VDS  (volts) 

Fig.  5.2.  MOSFFT  pundit lirough  characteristics  using  PISCTS,  GFM1X1, 
and  ('ADOPT.  The  current  path  in  the  .4 /tin  junction  device  is  at 
the  semiconductor  surface  for  biases  below  7V  and  in  the  bulk  for 
higher  biases.  The  . -la //in  junc!  ion  device  results  are  from  PISCThS 
only. 


crease  in  punehthrough  current  with  a  \2°c  increase  in  junction  depth  is  a 
rather  startling  result  and  illustrates  the  utility  of  numerical  simulation  in 
device  design.  'Phis  extreme  sensitivity  and  similar  sensitivities  to  junction 
curvature  and  impurity  gradients  make  anlvtica!  modeling  of  these  effects 
extremely  didicult,. 

The  three  simulation  programs  used  to  analyze  the  .1 //m  junction  depth 
device  agree  quite  well  throughout  the  range  of  drain  biases.  The  reason 
for  the  discrepancy  in  the  CFMINI  solution  at  Vps  —  61'  is  unknown  but 
may  be  related  to  the  shifting  of  the  dominant  punehthrough  current  path 
from  the  surface  to  the  bulk.  The  CADDFT  program  shows  a  variation  from 
the  other  two  at  drain  voltages  above  8VT  The  largest  variation  between  the 


throe  programs  is  approximately  a  factor  of  two.  This  is  larger  than  one 
would  like  but  certainly  sufficient  for  the  accuracies  typically  achieved  for 
subthreshold  device  characteristics.  Similar  comparisons  between  PISCES 
and  CADDET  for  devices  operating  above  threshold  show  excellent  agreement 
with  variations  typically  less  than  ;V  /.  Note  that  GEMINI  is  applicable  only 
in  the  subtheshold  and  low-drain-bias  linear  regions  of  operation. 

Inevitably,  differences  in  results  from  the  three  programs  can  be  traced  to 
differences  in  grid  placement,  unrealistic  physical  assumptions  for  the  simu¬ 
lated  conditions,  or  loosely  converged  solutions.  The  automatic,  grid  genera¬ 
tion  of  CADDKT  and  the  inflexibility  of  the  rectangular  grid  occasionally 
result  in  poorly  placed  grid  points.  The  limited  number  of  allowed  grid  points 
in  PISCES  also  may  result  in  sub-optimum  grid  placement.  The  CADDET 
assumption  of  a  one-dimensional  electric  field  in  the  oxide  and  the  GEMINI 
assumption  of  constant  quasi-Fermi  levels  can  cause  difficulty  if  the  device 
structure  or  bias  exceed  valid  ranges.  Under-converged  solutions  are  par¬ 
ticularly  noticeable  in  CADDET  due  to  the  solution  methods  and  convergence 
criteria  used. 

All  of  these  programs  arc  more  limited  in  accuracy  by  their  models  of 
higher  order  physical  phenomena,  however,  than  by  their  numerical  methods. 
Accurate  simulators  must  include  models  of  phenomena  such  as  velocity 
saturation,  field  and  concentration  dependent  mobility,  handgap  narrow¬ 
ing,  degenerate  statistics,  surface  states,  two-dimensional  impurity  profiles, 
Schoflky  barriers,  etc.  The  limitation  is  nol  one  of  implementation,  but 
rather  one  of  obtaining  an  accurate  model.  Much  current  controversy  sur¬ 
rounds  such  topics  as  handgap  narrowing  and  mobility  models.  Tin'  utility  of 
numerical  simulation  in  investigating  such  models  is  displayed  in  the  following 
section. 
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5.2  Strong  Inversion  Mobility 


The  subject  of  channel  mobility  under  strong  inversion  has  had  con¬ 
siderable  attention  in  recent  years  [5.3  5.5].  At  issue  is  the  variation  of  sur¬ 
face  mobility  with  substrate  doping,  vertical  electric  field,  crystal  orientation, 
substrate  bias,  and  interface  charge. 

Device  designers  have  traditionally  used  empirical  relations  for  predict¬ 
ing  device  performance  based  on  the  observ'd  effects  of  substrate  bias  and 
impurity  level  on  channel  mobility.  It  was  generally  understood  that  some 
form  of  scattering  caused  the  mobility  variations  but  that  the  effects  could  be 
parameterized  in  terms  of  the  substrate  bias  or  impurity  levels.  More  recently, 
it  has  become  evident  that  the  scattering  occurs  at  the  semiconductor- 
insulator  interface  and  that  the  shape  of  the  inversion  layer  charge  profile  is 
correlated  with  the  mobility  variations,  he.  the  mobility  reduction  is  greatest 
when  the  centroid  of  the  inversion  layer  profile  is  nearest  the  surface.  This 
compression  of  the  inversion  layer  at  the  surface  is  related  to  the  magnitude 
of  the  surface  electric  field,  thus  empirical  models  using  the  surface  field  as  a 
parameter  have  emerged  [5 .(>] .  The  observed  variations  with  subst  rate  doping 
and  bias  are  caused  by  the  difference  in  surface  fields  required  for  equivalent 
inversion  levels  at  the  different  substrate  impurity  levels  and  biases.  Even 
better  agreement  between  model  and  measurement  over  a  wide  range  of  sub¬ 
strate  conditions  has  been  achieved  with  use  of  the  “average”  or  effective  field 
in  the-  inversion  layer  [5.5]. 

From  Gauss’  law,  the  surface  electric  field  is  given  by 

r  Qi  b  Qn 

£  3 

while  the  effective  field  is  given  by 


(5.2) 


Eeff  = 


Q//2  +  Qu 

£3 


where  Qi  is  the  inversion  layer  charger  per  unit  area,  Qn  is  the  hulk  depletion 
charge  per  unit  area,  and  ta  is  the  semiconductor  permittivity.  The  current 
belief  is  that  phonon  and  surface  roughness  scattering  mechanisms  dominate 
the  room  temperature  surface  mobility  in  strong  inversion  and  phonon  and 
coulombic  (interface  and  oxide  charge)  scattering  in  weak  inversion;  however, 
the  quantification  of  these  effects  is  not  well  understood. 

The  point  of  the  foregoing  discussion  is  that  the  empirical  and  analyti¬ 
cal  models  for  surface  mobility  used  by  device  designers  and  in  circuit 
simulation  programs  necessarily  avoid  handling  the  basic  underlying  physi¬ 
cal  mechanisms  in  order  to  obtain  computationally  manageable  models. 
Numerical  simulation,  on  the  other  hand,  can  more  easily  accomodate  the 
underlying  physics  and  thus  can  be  used  by  device  designers  for  more  ac¬ 
curate  and  operating-region-independent  results,  and  by  device  physic:-  s  for 
studying  the  validity  of  their  models.  In  this  section,  PISCES  is  used  to 
evaluate  a  surface  mobility  model  in  which  the  mobility  varies  with  distance 
from  the  surface. 

A  very  long  channel  device  was  chosen  for  this  study  in  order  to  minimize 
the  influence  of  the  source  and  drain  on  the  channel  potential.  The  important 
device  parameters  are: 


oxide  thickness 

1000  A 

channel  length 

50  /nn 

substrate  doping 

1.2  X  10' 

fixed  interface  charge 

101 1  cm- 

It  is  well  known  that  the  carrier  mobility  varies  from  its  value  in  the 
bulk  to  a  lower  value  at  or  near  the  surface;  however,  the  value  of  the  surface 
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mobility  and  the  functional  form  of  the  variation  with  distance  are  not  known. 
A.  Gniidinger  and  II.  Tally  [5.7]  have  computed  the  thickness  of  inversion 
layers  to  be  on  the  order  of  100A  and  have  shown  quantum  mechanically  that 
the  carrier  density  peaks  at  distances  on  the  order  of  25  A  from  the  surface. 
It  is  reasonable  to  assume  that  the  variation  with  distance  is  monotonic,  thus 
the  mobility  model  chosen  is  one  in  which  the  mobility  decreases  exponentially 
with  distance  from  the  surface  with  characteristic  length  on  the  order  of  25 A. 
Mathematically,  it  is  expressed  as 

Ky)  —  Hb  -  {fib  -  H»)e~y/a  (5.3) 

where  is  the  bulk  mobility,  fis  is  the  mobility  exactly  at  the  surface,  y 
is  the  distance  from  the  surface,  and  o  is  the  characteristic  length  of  the 
variation.  Figure  5.3  shows  the  results  of  simulations  with  three  different 
values  of  the  model  parameters.  In  (a)  there  is  no  mobility  variation  with 
depth  and  the  classical  straight  line  variation  of  I D  with  V(;  is  observed.  In  (b) 
the  mobility  is  1000  cm2/V-s  in  the  hulk,  100  cm2/V-s  at  the  surface  and  has 
a  characteristic  length  of  50  A.  In  (c)  the  bulk  mobility  remains  the  same  but 
the  surface  mobility  is  reduced  to  10  cm2/V-s  and  the  characteristic  length 
to  33  A.  In  order  to  accurately  quantize  these  variations,  the  grid  spacing 
perpendicular  to  the  surface  is  very  small,  starting  at  less  than  10  A.  This 
small  spacing  is  somewhat  restrictive  since  it  consumes  large  numbers  of  grid 
points  and  could  cause  numerical  errors  in  the  difference  equations. 

The  reduced  drain  current  for  cases  (b)  and  (c)  of  Figure  5.3  reflects 
the  reduced  channel  mobility.  The  flattening  out  of  the  curves  at  higher 
gate  biases  is  characteristic  of  ICFFT’s  and  results  from  the  crowding  of  the 
inversion  layer  charge  closer  to  the  surface.  Figure  5.1  shows  the  relation 
between  the  Ip/Va  curves,  the  field  effect  mobility  (p7, ■/.;),  and  effective 
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Fig.  5.3.  Rolloff  of  drain  current  with  gate  bias  showing  effect  of  distancc- 
from-thc-suracc  rnobilily  model.  The  mobility  is  in  cm2/V-s  and 
y  is  in  A. 


mobility  (/<<•//).  These  are  expressed  mathematically  as 
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and 
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(5.5) 


where  L  and  W  arc  the  channel  length  and  width,  y(t  and  r;m  are  the  drain 
conductance  and  trarisconductanco,  qNinv  is  the  total  induced  charge  in  the 
channel  per  unit  area,  and  Cq  is  the  gate  capacitance  per  unit  area.  Note 
from  the  figure  that  these  two  mobilities  are  equal  at  the  steepest  part  of  the 
curve.  These  relations  arc  described  by  S.C.  Sun  [5.8]  whose  measurements 
are  used  for  t  he  comparisons  in  the  remainder  of  this  section. 


10G 


Fig.  5.4.  Definition  and  relationship  between  effective  mobility  fi,j/  and  field 
effect  mobility  fipE- 

Figure  5.5  shows  the  comparison  of  simulation  versus  measured  results 
for  the  device  described  earlier.  The  mobility  model  parameters  used  are 
fib  =  1286  cm2/V-s,  /xs  =  400  cm2/V-s,  and  a  —  50  A.  The  bulk  mobility 
value  is  chosen  based  on  the  substrate  impurity  concentration.  The  solid 
and  dashed  lines  represent  the  simulated  data  for  fi,.j /  and  fipv  respectively 
and  the  circles  and  squares  the  measured  values.  Although  the  simulated 
and  measured  results  arc  olTset,  they  have  the  same  shape.  Obviously,  the 
distance- from-the-surfacc  mobility  model  has  the  proper  effect  on  the  effective 
mobility  although  the  functional  form  or  parameter  values  of  the  model  may 
not  be  correct.  It  would  be  of  interest  to  fit  the  two  unknown  parameters,  //, 
and  a,  to  additional  measured  results  in  order  to  determine  the  parameter 
sensitivities  to  device  fabrication  and  structure  variations.  The  curve  fitting 
would  be  an  empirical  study;  however,  the  parameters  have  a  physical  basis 
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Fig.  5.5.  Comparison  of  simulated  and  measured  values  of  //,//  and  ///.•/.; 
with  varying  gate  voltage. 


and  could  shed  light  on  mobility  phenomena  as  opposed  to  analytical  models 
such  as  those  used  in  circuit  simulation  program  device  models  which  have 
no  physical  basis.  This  study  is  not  attempted  since  it  is  not  the  intent  of 
this  work  to  develop  a  new  mobility  model,  but  it  is  suggested  as  an  area  of 
further  research. 

The  maximum  value  of  mobility  occurs  at  the  point  labelled  in 

Figure  5.5.  Since  this  peak  always  occurs  at  a  gate  voltage'  just  above 
threshold,  the  surface  is  only  weakly  inverted.  As  substrate'  doping  increases, 
the  surface  liedel  required  for  the  same  degree  of  surface  inversion  also  in¬ 
creases  resulting  in  more  crowding  of  the  inversion  layer  charge  anel  thus  lower 
thnnx-  This  reduction  in  maximum  mobility  with  increased  substrate*  doping 
is  seen  in  Figure  5.6.  Both  measured  anel  simulated  results  are*  shown  along 
with  the  bulk  mobility  values  for  comparison.  The>  simulated  results  were  ed>- 
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Fig.  5.(i.  Comparison  of  simulated  and  measured  values  of  /imas  sensitivity 
to  substrate  doping. 

tained  at  substrate  dopings  of  1.2  X  1 0 1 5 ,  1.2  X  10H),  and  1.2  X  I017  cm-3. 
The  simulated  results  are  comparable  to  the  measured  values  but  show  a 
steeper  slope  indicating  that  the  surface  mobility  reduction  effect  is  too  strong 
with  the  chosen  parameters. 

One  of  the  principal  points  of  this  section  is  that  numerical  simulation 
can  be  useful  in  the  evaluation  of  physical  models.  Attempting  to  match  the 
measured  data  shown  here  and  available  elsewhere  by  modifying  the  mobility 
model  should  provide  additional  insight  into  the  characteristics  of  surface 
inversion  layer  mobility.  The  other  point  is  that  the  use  of  the  most  fun¬ 
damental  physical  models  possible  results  in  the  most  powerful  and  versatile 
device  simulation  programs.  The  relatively  simple  distance-from-the  surface 
mobility  model,  for  example,  is  easier  to  implement  than  an  effective  field 


mobility  reduction  model  and  much  more  universal  than  substrate  bias  and 
substrate  doping  mobility  reduction  models. 


5.3  Summary 


Two  examples  of  application  of  the  PISCES  program  are  presented- 
punchthrough  current  simulations  on  an  implanted  channel  NMOS  transistor 
and  evaluation  of  a  distance-from-the-surface  mobility  model.  The  PISCES, 
GEMINI,  and  CADDET  programs  are  compared  for  the  punchthrough 
simulations.  Their  results  agree  throughout  the  range  of  simulated  biases. 
Simulations  also  show  that  a  12%  increase  in  source/drain  junction  depth 
can  result  in  a  50  fold  increase  in  punchthrough  current.  Simulations  using 
the  depth  dependent  mobility  mode!  indicate  that  such  a  model  may  be  useful 
for  device  simulation  in  lieu  of  field  dependent  models.  An  exact  form  for  the 
model  is  not  pursued  but  is  suggested  as  an  area  of  further  research.  The 
effects  of  mobility  reduction  with  increased  gate  bias  and  increased  substrate 
doping  arc  demonstrated.  One  disadvantage  of  tne  depth  dependent  mobility 
model  is  the  need  for  very  fine  grid  spacing  normal  to  the  surface. 

The  next  chapter  summarizes  the  conclusions  and  recommendations  of 
this  work. 
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Chapter  6 


CONCLUSION 

The  two-dimensional  structure  of  modern  semiconductor  devices  demands 
the  use  of  two-dimensional  numerical  simulation  in  device  design.  The  ap¬ 
plication  of  analytical  models  simply  cannot  accurately  account  for  the  highly 
two-dimensional  impurity  profiles  and  potentials  and  their  interaction. 

All  of  the  two-dimensional  device  simulation  programs  currently  avail¬ 
able  are  too  restrictive  to  be  of  great  service  to  device  designers  as  witnessed 
by  the  slow  acceptance  of  the  few  commercially  available  device  simulation 
programs.  Restrictions  in  allowable  device  structures,  grid,  computation 
time,  memory  requirements,  accuracy  of  physical  models,  and  ease  of  user 
interface  are  seen  to  some  degree  in  all  programs.  Various  aspects  of  those 
limitations  have  been  addressed  in  this  work  through  the  development  and 
application  of  a  device  simulation  program,  PISCES.  The  program  contains 
some  of  the  same  limitations  but  is  extremely  flexible  and  allows  invest  igation 
of  many  of  the  device  simulation  program  restrictions. 


6.1  Summary 

The  allocation  of  grid  in  various  regions  of  a  device  has  been  addressed 
in  terms  of  grid  type,  grid  density,  device  structures,  and  impurity  profiles. 
These  analyses  show  that  the  densest  grid  is  required  in  regions  of  high  net 
charge  density,  large  gradients  of  net  charge  density  or  large  gradients  of 
potential.  The  use  of  reflecting  boundary  conditions  along  the  sides  of  the 

111 


device  is  shown  to  require  significant  lateral  extensions  of  the  source  and  drain 
in  order  to  accurately  represent  the  device  potentials.  Simulations  have  also 
demonstrated  the  inadequacy  of  rectangular  uniformly  doped  approximations 
to  the  source/drain  regions.  Even  a  very  coarse  approximation  to  a  Gaussian 
source/drain  profile  is  shown  to  provide  very  good  results. 

The  application  of  finite  difference  discretization  of  Poisson’s  equation 
and  the  current  continuity  equation  to  an  irregular  triangular  grid  has  been 
presented  including  the  special  cases  of  obtuse  triangles.  A  more  consistent 
area  allocation  scheme  has  been  presented  along  with  a  simple  technique  for 
avoiding  negative  coupling  coefficients  in  the  Poisson  discretization  for  obtuse 
triangles.  The  quasi-two-dimensional  discretization  of  the  continuity  equation 
using  the  Scharfotter-Gumrnel  algorithm  is  accompanied  by  a  proof  of  the 
non-existence  of  a  fully  two-dimensional  form. 

A  variety  of  matrix  solution  techniques  have  been  compared  in  terms 
of  their  applicability  to  device  simulation.  The  flexibility  and  stability  of 
the  LIJ  decomposition  method  are  offset  by  the  rapid  growth  of  solution 
time  and  memory  requirements  with  grid  size.  The  opposite  can  be  said  for 
the  iterative  methods  of  SfP  and  fCCG  which  do  not  grow  as  rapidly  with 
grid  size  but  are  more  restrictive  in  grid  type  and  sensitive  to  the  matrix 
coefficients  in  solution  time. 

Evaluation  of  the  LU  decomposition  method  has  shown  that  proper 
numbering  of  the  grid  can  result  in  time  and  memory  savings.  Numbering 
a  rectangularly  connected  grid  in  the  shortest  direction  (i.c.  row  or  column) 
minimizes  the  generation  of  non-zero  matrix  elements  and  thus  reduces  the 
matrix  equation  solution  time.  A  more  exotic  renumbering  scheme,  nested 
dissection,  is  shown  to  reduce  the  solution  time  from  0(i'V2)  to  0(/V'<//")  and 
storage  requirements  from  0(NA/2)  to  ()(/V  Inn).  These  results,  obtained  on 
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triangular  grids,  match  those  reported  for  rectangular  grids.  A  drawback  of 
the  nested  dissection  method  is  that  vector  computers  cannot  be  used  to  full 
advantage  because  of  the  shorter  vectors  generated  by  this  numbering. 

The  convergence  properties  of  Cumrnel’s  alternating  method  for  the  solu¬ 
tion  of  Poisson’s  equation  and  the  continuity  equation  has  been  thoroughly 
analyzed.  The  method  works  exceedingly  well  for  devices  biased  below 
threshold,  generally  requiring  less  than  five  iterations.  Simulations  of  devices 
biased  above  threshold  but  with  very  little  current  flowing  also  converge 
rapidly.  Significant  convergence  problems  occur,  however,  for  devices  biased 
in  the  linear  and  saturation  regions  of  operation  where  substantial  currents 
are  flowing. 

Four  different  methods  have  been  derived  in  the  course  of  this  work  for 
accelerating  the  convergence  of  the  alternating  method  in  these  operating 
regions-  projection  of  the  initial  guess  from  previous  solutions,  use  of  only 
one  Poisson  iteration  per  alternating  iteration,  overrelaxation  of  the  electros¬ 
tatic  potential  solutions  and  reduction  of  the  Poisson  linearizing  term.  The 
use  of  combinations  of  these  techniques  reduces  the  average  solution  time  by 
approximately  a  factor  of  four.  Using  these  acceleration  techniques,  as  the 
drain  bias  increases  the  solution  time  for  MOSPET  drain  bias  stops  is  seen  to 
increase  until  reaching  saturation  and  then  decrease  as  the  device  is  biased 
deeper  into  saturation.  This  solution  lime  reduction  with  increasing  bias  in 
the  saturation  region  is  seen  only  when  using  these  acceleration  techniques 
and  tends  to  counter  claims  that  the  alternating  method  is  not  practical  for 
simulation  of  devices  above  threshold. 

Even  with  the  acceleration  techniques,  the  convergence  above  threshold 
is  still  relatively  slow.  Plots  of  surface  potential  error  versus  location  along 
the  channel  show  a  single  sinusoidal  variation  of  the  error  from  source  to 
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drain.  The  sinusoid  undergoes  very  slow  decaying  oscillations  in  a  wave-like 
manner  as  the  iterations  progress.  This  first  harmonic  in  spatial  frequency 
of  the  surface  potential  error  appears  to  be  a  dominant  factor  in  the  slow 
convergence  of  the  alternating  method  above  threshold. 

Comparisons  of  the  PISCES,  GEMJAI  and  CAPDET  device  simulation 
programs  has  shown  agreement  to  within  a  factor  of  two  in  subthreshold 
simulations  and  to  within  a  few  percent  above  threshold.  These  and  oilier 
differences  between  simulation  program  solutions  can  usually  be  traced  to 
difference  in  grid  placement,  physical  assumptions  or  convergence  tolerance. 

A  depth  dependent  mobility  model  lias  been  implemented  in  order  to 
examine  its  feasibility  and  to  demonstrate  the  ability  of  numerical  simulation 
programs  to  use  more  fundamental  physical  models  than  the  analytical  or 
empirical  models.  Simulations  using  this  model  have  shown  reasonable  agree¬ 
ment  with  measurement  for  mobility  variations  with  gate  bias  and  substrate 
doping.  This  application  of  device  simulation  also  demonstrates  the  utility 
of  such  programs  in  the  evaluation  of  physical  models. 

6.2  Rccommt'ndntions 

As  with  most  scientific  endeavors,  there  appear  to  he  more  questions  at 
the  end  than  there  were  at  the  start.  The  PISCES  program  is  extremely 
well  suited  to  the  studies  performed  in  this  work  hut  is  unsatisfactory  as  a 
device  designers  tool.  Several  modifications  to  the  program  are  suggested 
witli  varying  degrees  of  additional  research  required. 

'Pile  rectangularly  connected  triangular  grid  of  PISCES  appears  to  be  a 
reasonable  compromise  in  terms  of  llexihiHity  in  matching  device  structures, 


ease  of  generation,  and  compatibility  with  matrix  equation  solution  tech¬ 
niques.  An  automatic  grid  generation  scheme  should  be  developed  for  im¬ 
plementing  this  grid  based  on  the.  grid  density  criteria  described  in  this  work. 
Grid  refinement  with  changing  bias  should  be  studied  with  a  consideration 
for  the  tradeoff  between  the  increased  accuracy  achieved  and  the  possibility 
of  slight  discontinuities  in  the  device  characteristics  due  to  the  changing  grid. 

The  LU  decomposition  matrix  solution  method  used  in  PISCPS  is  not 
the  optimum  method  for  device  simulation  and  should  be  replaced.  Studies 
should  be  performed  on  comparisons  of  the  SIP  and  1CGG  methods  as  they 
appear  to  be  the  best  suited  replacements.  Implementation  of  either  of  these 
iterative  matrix  solution  techniques  should  be  accompanied  by  study  of  the 
effectiveness  of  fully  merging  the  alternating  solution  iterations  by  mixing 
the  Poisson  and  continuity  equation  matrix  iterations. 

Although  the  obtuse  triangle  discretization  is  relatively  sound  and  many 
apparently  successful  simulations  have  been  performed  with  obtuse'  triangles 
in  the  grid,  their  effect  is  not  fully  understood.  Some  high  resolution  contour 
plots  in  regions  with  large  numbers  of  highly  obtuse  triangles,  for  example, 
have  shown  slight  distortions  in  the  cquipotential  contours.  Additional  re¬ 
search  is  suggested  to  quantify  these  effects. 

The  slow  convergence  of  the  alternating  method  for  simulation  of  devices 
above  threshold  remains  a  problem  in  spile  of  the  acceleration  achieved  in 
this  work.  The  significance  of  the  oscillations  seen  in  the  surface  potential 
should  be  investigated  with  the  possibility  of  supplementing  the  Poisson  and 
continuity  equation  coupling  along  the  length  of  the  channel. 

A  more  definitive  study  of  the  depth  dependent  mobility  model  should  be 
performed  in  order  to  determine  the  proper  functional  form  and  parameter 
values  for  the  model  and  its  practicality  for  use  in  device  simulation.  The 
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grid  density  requirements  for  implementation  of  the  model  should  also  he 
considered.  Other  poorly  understood  phenomena  such  as  weak  inversion 
mobility  should  also  be  examined  via  device  simulation. 


Appendix  A 


PISCES  DEMONSTRATION  EXAMPLE 

Parts  of  the  PISCES  program  have  been  described  in  the  text  of  this 
work;  however,  sufficient  detail  for  a  thorough  understanding  of  the  program 
functions  is  not  provided.  This  appendix  provides  the  necessary  detail  in 
the  form  of  an  example  which  demonstrates  most  of  the  program  features. 
The  device  simulated  is  an  N-channel  Silicon  MOS1ET  with  the  following 
parameters: 


gate  length 

1 .5  pm 

gate  oxide 

500  A 

field  oxide 

4000  A 

junction  depth 

.4  pm 

substrate  doping 

2  X  10' 5  cm-3 

channel  implant  dose 

3.1  X  10"  cm-2 

fixed  interface  charge 

1  X  10'°  cm.-'2 

gate  material 

n-type  polysilicon 

source/drain  contacts 

aluminum. 

Figures  A.  I  and  A. 2  show  the  input  deck  for  the  program.  The  deck  is 
divided  into  two  parts  to  demonstrate  the  saving  and  restoring  of  program 
data  files.  Csually,  one  would  split,  the  deck  into  several  pieces  so  that  each 
step  in  the  simulation  may  be  verified  before  proceeding  to  the  next. 

The  first  item  in  each  line  is  the  card  name  and  the  remaining  items  are 
parameters.  There  are  three  types  of  parameters  numeric,  alphanumeric, 
and  logical.  Numeric  parameters  are  followed  by  an  equal  sign  and  a  numeric 
value.  Alphanumeric  parameters  are  also  followed  by  an  equal  sign  but.  may 


TITLE  MOSFET  EXAMPLE 

$  ***  Generate  mesh  *** 

MESH  RECTANGULAR  NX-31  NY-22  OUTF I LE -MESH1  BULKDOP-2E15  P.TYPE 

X . MESH  NODE  =  1  LOCATION-1  RATIO-1 

X.M  N-3  L  =  1 . 5  R  = . 7 1 

X.M  N-9  L  =  1 . 9  R=  .  8 

X.M  N=  15  L=2 . 25  R-l. 11 

X.M  N-2 1  L-2.6  R-.9 

X.M  N-27  L-3  R-l. 25 

X.M  N  =  29  L  =  3 . 5  R-l. 414 

X. M  N-3 1  L-4.5  R-l. 414 

Y. M  N  =  1  L  =  - .05  R-l 

Y.M  N-4  L  =  0  R-l 

Y.M  N-22  L=3  R-l. 25 

$  ***  Expand  field  oxide  *** 

SPREAD  LEFT  WIDTH-. 5  UPPER-1  LOWER-4  THICKNESS-. 4  ENCROACH-1 
+  VOL. RAT-. 4 

SPR  RIGHT  W-1.5  UP-1  LO-4  THICK-. 4  ENCR-1  V0L-.4 

$  *  *  *  Match  junctions  *  *  * 

SPR  LEFT  W-.8  UP-4  L0-10  Y-.41  ENCR-.9  GRADING-. 7 

SPR  RIGH  W-1.8  UP-4  L0-10  Y-,41  EN-.9  GR-.7 

$  ***  Identify  insulator  and  semiconductor  regions  *** 

REGION  NUMBER-1  X. LOW-l  X. HIGH-31  Y. LOW-1  Y. HIGH-4  INSULATOR 
REG  NUM-2  X.L-1  X.H-31  Y.L-4  Y.H-22  SEMICONDUCTOR 

$  ***  Dope  the  semiconductor  *** 

$  ***  substrate  *** 

DOPING  P.TYPE  CONCENTR-2E 15  UNIFORM 
$  ***  channel  implant  *** 

DOP  P  DOSE-3. 4E11  Y. PEAK-. 18  Y . CHARAC- . 2404  GAUSSIAN 

$  ***  source  and  drain  *** 

DOP  DONOR  C0NC-4E  19  LEFT . JUN  Y.JUNC-.4  Y. PEAK-0  GAUSS 

+  X. RIGHT-1. 5  XY. RATIO-1 

DOP  DONOR  CONC-4E19  RIGHT. J  Y.J-.4  Y.P-0  GAUSS  X.L-3  XY-1 

$  ***  Fixed  surface  states  *** 

QF  CONCENTR  =  IE  10  X. LOW-3  X. HIGH-27  Y. LOW-4  Y. HIGH-4 

$  ***  Identify  electrode  locations  *** 

ELECTROD  NUMBER-1  X. LOW-3  X. HIGH-27  Y. LOW-1  Y.HIGH-l 
ELEC  N  =  2  X.L-1  X.H-31  Y.I.-22  Y.H-22 

ELEC  N-3  X.L-1  X.H-2  Y.L-4  Y.H-4 
ELEC  N-4  X.L-28  X.H-31  Y.L-4  Y.H-4 

$  ***  Print  vertical  grid  info  *** 

PRINT  POINTS  IX. MIN-15  IX. MAX-15 

$  **»  Plot  grid  and  junctions  *** 

FLOT.2D  X. MIN-1  X. MAX-4. 5  Y.MIN--.3  Y. MAX-3  NO. TOP  BOUNDARY 
+  JUNCTION  GRID 

$  ***  End  *** 


Fig.  A.  1 .  Sample  PISC'FS  input,  card  dock  for  mesh  generation  and  device 
structure  definition. 
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TITLE  PERFORM  SOLUTIONS 
S  ***  Get  mesh  *** 


MESH  IN=MESH1 

$  ***  Perform  symbolic  matrix  factorization  *** 
SYMB  OUT  =  SYMB1 


$  ***  Prepare  for  initial  solution  *** 

SETUP 

IN  IT  PRINT  TEMPERAT  =  300  P.ELECT  =  2 

$  • 

specify  materials  *** 

MATERIAL 

NUMBER1 1  OXIDE 

MATER 

NUM=  2  SILICON 

$  •** 

specify  contacts  *** 

CONTAC 

NUMBERS  N  .  POLY 

CONTAC 

NUM=2  NEUTRAL 

CONTAC 

NUM=3  ALUMINUM 

CONTAC 

NUM=4  ALUM 

$  *** 

Specify  mobility  models  *** 

MOBILITY 

VSAT  C0NM0B 

S  ***  Solve  initial  solution  *** 

SOLV 

PRINT  OUT  =  EXOUT 0 

$  ***  Step 

i  gate  bias  *  *  * 

SETUP 

INF  =  EXOUT 0  PREVIOUS  V1=0 

SOLVE 

PRINT 

SETUP 

PROJECT  Vl=2 

SOLVE 

0UTFILE=EX0UT1  PRINT 

$  ***  Step  drain  bias  *** 

SETUP  PREVIOUS  V4=  .  5 

SOLVE  SINGLE  ACCEL  OUT=EXOUT2  PRINT 

SETUP  PROJECT  VSTEP= . 5  NSTEPS=3  ELEC=4 
SOLVE  SINGLE  ACCEL  OUT=EXOUT3  PRINT 

$  ***  Plot  results  at  Vg  =  2,  Vd  =  2  *** 

PLOT. 2  X.MIN=1  X . MAX=4 . 5  Y.MIN=-.3  Y.MAX=3  BOUND  NO. TOP  JUNC  DEPL 
CONTOUR  POTENTIAL  MIN.VAL=-.2  MAX.VAL=i.6  DEL.VAL=.2 
PLOT. 2  X . MIN= 1  X . MAX=4 . 5  Y.MIN=-.3  Y.MAX=3  BOUND  NO. TOP  JUNC  DEPL 
CONTO  QF.POT  MIN=.2  MAX=2  DEL=.2 

$  ***  End  *** 


Fig.  A. 2.  Sample  PISCES  input  card  deck  for  specifying  device  material 
characteristics  ami  obtaining  simulation  solutions. 


have  any  alphanumeric  character  as  a  value.  Logical  parameters  may  he 
followed  by  an  equal  sign  and  the  words  true  or  false  or  may  appear  alone  in 
which  case  they  are  assigned  a  logical  value  of  true.  A  in  the  first  column 
indicates  a  continuation  of  the  previous  line.  Note  that  either  card  names 
or  parameters  names  may  be  shortened  if  the  resulting  name  is  unambigious. 
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The  card  and  parameter  names  recognized  by  the  program  (8  characters 
maximum)  are  shown  in  upper  case  letters  with  the  remainder  of  the  name  in 
lower  case  letters  for  clarity.  The  remainder  of  this  appendix  details  the  use 
of  cards  and  parameters  by  describing  their  use  in  the  sample  input  decks  of 
Figures  A.l  and  A. 2. 

TITLE 

The  TITLE  card  has  no  parameters.  All  of  the  characters  after  the  card 
name  are  stored  and  used  as  a  header  for  all  printed  listings. 

$  or  COMMENT 

Either  $  or  COMMENT  may  be  used  to  specify  a  comment  line  which 
is  ignored  by  the  program. 

MESII 

The  MESH  card  indicates  the  beginning  of  a  sequence  of  cards  serving 
to  define  the  device  structure.  The  sequence  is  terminated  when  a  non-mesh- 
defining  card  is  encountered.  Most  of  the  cards  must  appear  in  the  order 
given  in  order  to  properly  define  the  device.  A  KECTANGUlar  mesh  (grid)  is 
specified  meaning  that  the  grid  nodes  will  initially  lie  at  the  intersections  of 
parallel  horizontal  and  vertical  lines.  Distortion  of  this  grid  is  allowed  later. 
There  are  31  vertical  grid  lines  (NX)  and  22  horizontal  grid  lines  (NY).  At 
the  termination  of  the  mesh  sequence,  all  of  t  he  structure  data  will  be  st  ored 
in  a  file  (OUTFILE)  called  MESII  1.  The  substrate  doping  (BULKDOP)  is 
2  X  1015  cm-3  and  is  p-type  (P.TYl’E).  The  substrate  may  also  be  specified 
as  N.TYPE.  If  a  mesh  file  has  been  previously  stored,  all  structure  data  may 
be  read  with  the  single  parameter  INFILE  and  the  name  of  the  iile. 
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XJvJESIJ 

The  X.MESH  card  specifies  the  location  along  the  x  axis  of  one  of  the 
vertical  grid  lines.  The  first  NODE  (actually  all  nodes  on  the  first  vertical 
line)  has  an  x-axis  coordinate  (LOCATION)  of  1  fun.  The  origin  of  the  x  axis 
may  be  arbitrarily  chosen.  The  RATIO  value  has  no  meaning  for  the  first 
node. 

The  next  card  is  also  an  X.MES1I  card  but  the  card  name  and  parameter 
names  are  shortened  for  ease  of  typing.  All  of  the  cards  in  this  example  follow 
this  same  pattern  in  which  full  card  and  parameter  names  are  used  on  the 
first  occurrence  of  a  card  type,  but  shortened  names  are  used  thereafter. 
The  s  ond  X.MES1I  card  places  node  3  at  1.5  //m.  The  RATIO  parameter 
specific  that  the  spacing  between  vertical  grid  lines  2  and  3  should  be  only 
.71  of  the  spacing  between  lines  1  and  2.  If  there  were  more  grid  lines  specified 
in  the  interval,  then  each  successive  space  (from  left  to  right)  would  be  .71 
as  large  as  the  previous  space.  The  additional  X.MESI1  cards  specify  the; 
remaining  grid  lines,  locations  and  spacing  ratios  up  to  the  rightmost  edge  of 
the  simulation  region  at  -1.5 /zm. 


Y.MES1I 


The  Y.MES1 !  card  serves  the  same  function  as  the1  X. MINI  I  card  but 
in  the  orthogonal  direction.  The  y  axis  is  positive  downward.  The  first 
horizontal  grid  line  is  placed  at  the  top  of  the  gate  oxide.  The  fourth  grid 
line  is  placed  at  the  oxide-semiconductor  interface  and  is  chosen  as  the  origin. 
The  last  grid  line  is  placed  at  the  bottom  of  the  simulation  region,  3 /cm  deep 
into  the  substrate. 
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SPREAD 

The  SPREAD  card  is  used  to  distort  the  grid  in  the  vertical  direction 
in  order  to  match  device  surface  or  interface  shapes  or  other  internal  device 
structure.  The  operation  results  in  horizontal  grid  lines  which  are  vertically 
displaced  on  either  the  left  or  right  side  of  the  device  with  a  smooth  variation 
of  this  displacement  across  the  device.  The  first  two  spread  cards  expand  the 
oxide  region  on  the  LEFT  and  RIGHT  sides  of  the  device  from  the  gate  oxide 
thickness  to  the  field  oxide  THICKNESs  of  4000  A.  The  second  two  distort 
the  grid  near  the  semiconductor  surface  to  match  the  source/drain  junction 
profiles.  The  first  spread  card  expands  the  grid  between  lines  1  (UPPER) 
and  4  (LOWER)  specified  later  as  the  oxide  region,  for  a  WIDTH  of  .5  /an 
from  the  LEFT  edge  of  the  device.  The  ENCROACHment  factor  specifies  the 
abruptness  of  the  transition  from  spread  to  non-spread  grid.  The  VOL.RATio 
parameter  specifies  the  ratio  of  the  downward  displacement  of  the  lower  grid 
line  to  the  net  increase  in  thickness,  corresponding  to  the  volume  ratio  of 
silicon  consumed  to  oxide  grown  in  thermal  oxidation. 

In  the  third  and  fourth  SPREAD  cards,  the  UPPER  grid  line,  the  oxide- 
semiconductor  interface,  is  not  moved  but  the  LOWER  line  is  moved  to 
the  Y.LOWEI?  coordinate  of  .41  /mi  which  is  just  below  the  source/drain 
junctions.  The  spacing  of  all  grid  lines  in  between  is  changed  to  a  GRADING 
of  .7  in  order  to  provide  t  he  proper  grid  placement,  on  the  steep  source/drain 
impurity  profiles.  The  GRADING  value  is  used  exactly  like  the  RATIO 
parameter  in  the  X.MESll  card. 

REGION 

The  REGION  card  is  used  to  define  the  INSULATOr  and  SEMK'ONDuc- 
tor  regions  of  the  device.  It  can  also  be  used  to  rigidly  restrict  the  region  of 
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the  device  which  receives  impurity  doping.  Each  region  must  he  sequentially 
numbered  using  the  NUMBER  parameter.  The  X.LOW,  X.IIIG1I,  Y.LOW, 
and  VMIIGH  parameters  specify  the  grid  lines  which  bound  the  region.  In 
this  typical  example,  the  entire  semiconductor  substrate  is  contained  in  one 
region,  and  the  oxide  in  another. 

DOPING 

The  DOPING  card  is  used  to  add  impurities  to  the  device  within  the 
bounds  of  the  most  recent  REGION  card.  The  first  DOPING  card  specifics 
the  substrate  doping  to  be  a  UNIFORM  distribution  of  P.T5  PE  impurity 
with  a  CONCENTRalion  of  2  X  10 15  cm"3. 

The  second  DOPING  card  specifies  the  channel  implant  as  a  GAUSSIAN 
implant  of  a  I’. TYPE  impurity  with  a  DOSE  of  3.4  X  10n  vm~2,  a  peak 
(Y.PEAK)  at  .18  pm  and  Cl  lARACTEristic  length  (y/'2<r)  of  .2101  pm. 

The  third  DOPING  card  specifies  a  GAUSSIAN  source  doping  with  a 
peak  CONGENTRation  of  4  X  1019  cm-3  DONOR  impurities,  with  the  peak 
(Y.PEAK)  at  the  origin  and  the  junction  (Y..JUNCT1)  at  .1  /mi  computed 
using  the  background  doping  on  the  left  side  (LEKT..JUN)  or  the  device.  The 
impurity  distribution  is  specified  to  be  uniform  in  the  lateral  direction  from 
the  left  edge  of  the  device  (by  default)  to  the  1.5  pm  location  on  the  x-axis 
(X. RIGHT).  This  point  corresponds  to  the  right  edge  of  a  diffusion  or  implant 
window.  The  lateral  profile  beyond  the  1.5pm  coordinate  is  also  Gaussian 
(by  default)  and  the  characteristic  length  in  the  x  direction  is  specified  bv 
XY. RATIO  to  be  equal  to  the  eharaef eristic  length  in  the  y  direction,  resulting- 
in  cylindrical  junctions. 

The  fourth  DOPING  card  specifies  the  same  doping  profile  for  the  drain. 
Additional  parameters  allow  the  impurity  type  to  be  specified  as  .N.'IA  PE  or 
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ACCEPTOR  and  permit  a  complementary  error  function  lateral  impurity 
profile  by  specifying  X.ERFC. 

QF 

The  QF  card  is  used  to  specify  the  fixed  surface  states  charge  at  an 
insulator-semiconductor  interface.  The  CONCENTRation  is  1  X  1010  crn-2 
and  exists  along  horizontal  grid  line  number  4  (Y.LOVV,  Y.HIGII)  from  ver¬ 
tical  grid  line  number  3  (X.LOW)  to  vertical  grid  line  number  27  (X.111GH). 

TCLECTRODe 

The  ELECTRODe  card  specifics  nodes  in  the  grid  at  which  the  potential 
boundary  conditions  will  be  applied.  Each  group  of  nodes  is  assigned  a 
NUMBER  which  is  used  to  reference  the  group.  The  node  clusters  are 
specified  by  the  bounding  grid  lines  as  in  the  QF  card  (X.LOW,  X.1I1GH, 
Y.LOW,  Y.HIGII).  Electrode  number  1  is  the  gate,  2  is  the  substrate,  3  is 
the  source  and  4  is  the  drain. 

This  is  the  last  mesh-defining  card.  Reading  the  next  card  terminates 
the  MESH  sequence  and  causes  a  final  computation  of  mesh  parameters  and 
stores  the  mesh  data  in  the  specified  output  file. 

PRINT 

The  PRINT  card  provides  terminal  or  fine  printer  listings  of  a  large 
variety  of  simulation  information.  The  information  is  printed  only  for 
areas  of  the  device  within  a  window.  The  window  may  be  specified  by 
providing  the  grid  line  boundaries  (IX. MIN,  IX. MAX,  JY.MIN,  IY.MAX)  or 
coordinate  boundaries  (X.MIN,  X.MAX,  Y.M1N,  Y.MAX).  The  POINTS 
parameter  prints  data  associated  with  nodes  in  the  grid.  Other  allow- 
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able  PRINT  parameters  are:  ELEMENTS — nodes  composing  each  triangle; 
GEOMETRY — coupling  coefficients;  SOLUTION — potential,  quasi-FERMI 
potential  and  electron  concentration  at  each  node;  MATERIAL — detail  on 
the  material  parameters  for  each  device  region. 

PLOT. 2D 

The  PLOT. 2D  card  makes  a  two-dimensional  plot  of  specified  device 
characteristics.  The  plot  window  is  specified  by  coordinates  (X.MIN,  X.MAX, 
Y.M1N,  Y.MAX),  which  may  lie  inside  or  outside  of  the  simulation  region. 
The  BOUNDARY  parameter  plots  the  device  external  boundary,  interfaces 
and  electrodes,  JUNCTION  plots  the  metallurgical  junctions,  GRID  plots  the 
triangular  grid,  and  DEPL.EDGe  plots  the  depletion  edges.  NO. TOP  inhibits 
the  plotting  of  tic  marks  across  the  top  of  the  plot  and  NO. TIC  inhibits  al! 
tic  marks.  NO. CLEAR  inhibits  clearing  of  the  display  between  plots  to  allow 
superimposed  plots.  A  plot  file  may  be  generated  by  specifying  OUTPILE 
and  a  file  name.  The  resulting  plot  is  shown  in  Figure  A. 3.  Further  plotting 
capability  is  provided  by  the  CONTOUR  card  to  be  described  later. 

Figure  A. 2  shows  the  card  sequence  which  performs  the  solutions  on  the 
device.  In  this  card  sequence,  the  MESII  card  merely  reads  the  file  MES1I  l 
which  contains  all  of  the  necessary  device  structure  information. 

.SYMBOLIC 

The  SYMBOLIC  card  invokes  the  symbolic  factorization  of  the  coefficient 
matrix  for  the  LU  decomposition.  INFILE  and  OUTITLE  are  used  to  read 
or  write  the  pointer  arrays  resulting  from  the  factorization.  The  DISSECT 
parameter  may  be  used  to  perform  a  nested  dissection  renumbering  on  the 


125 


Fig.  A. 3.  Plot  of  the  mesh  generated  by  the  PISCES  example. 

grid.  Alternatively,  MINIMIZE  may  be  used  to  renumber  the  grid  by  rows 
if  there  are  fewer  nodes  across  a  row  than  down  a  column.  The  PRINT 
parameter  prints  a  summary  of  relevant  factorization  parameters. 

SETUP 

The  SETUP  card  computes  the  coefficient  matrix  and  the  initial  guess 
prior  to  every  solution.  It  may  be  followed  by  a  sequence  of  cards  which 
specify  or  modify  various  device  parameters.  The  sequence  is  terminated 
when  a  non-setup-sequence  card  is  encountered.  Since  this  is  the  first  solution 
performed  on  this  device  the  INITIAL  parameter  is  specified.  This  resits  in 
a  charge- neutral  initial  guess  and  a  flat-band  bias  assignment.  The  PRINT 
parameter  invokes  a  listing  of  the  SETUP  parameter  values  on  termination  of 
the  SETUP  sequence.  The  device  TEMPERATure  is  specified  to  be  300 °K . 
The  P.ELECTRode  parameter  assigns  the  substrate  bins  (electrode  number 
2)  as  the  hole  quasi-Fermi  level  value.  An  alternative  way  to  set  the  hole 
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quasi-Fermi  level  is  to  explicitly  specify  it  with  the  P.BIAS  parameter. 

In  the  second  occurrence  of  the  SETUP  card  (several  cards  down  the  list), 
the  PREVIOUS  solution  is  used  as  the  initial  guess  for  the  next  solution.  The 
previous  solution  is  read  using  INFILE.  Actually,  reading  the  solution  from 
the  file  was  not  necessary  since  the  two  most  recent  solutions  are  always  stored 
in  memory.  INFILE  reads  a  stored  solution  into  the  most- recent- solution 
array  and  IN2FILE  reads  into  the  second-most- recent  solution  array.  The 
device  electrode  bias  levels  are  set  using  VI  through  V9  and  VTEN  where 
the  number  corresponds  to  the  electrode  number.  Here  the  gate  is  set  to 
zero  volts.  All  electrode  voltages  not  explicitly  set  are  kept  at  their  previous 
values. 

The  third  occurrence  of  the  SETUP  card  demonstrates  the  use  of  the 
PROJECTion  parameter  for  extrapolating  an  initial  guess  from  two  previous 
solutions.  Only  one  electrode  bias  is  allowed  to  change  between  the  new  and 
two  previous  solutions.  Here  the  gate  (electrode  1)  has  previous  values  of 
flat- band  and  zero  volts  and  a  new'  value  of  two  volts. 

The  fourth  SETUP  card  sets  the  drain  voltage  to  .5  volts. 

The  fifth  SETUP  card  demonstrates  the  bias  stepping  capability  of 
PISCES.  VSTEP  sets  the  bias  step  size,  NSTEPS  sets  the  number  of  bias 
steps,  and  ELECTRODe  specifies  the  electrode  number  (the  drain)  being 
varied.  When  bias  steps  are  specified  in  this  way,  only  one  SETl  ’P/SOLVE 
combination  is  required  for  the  range  of  bias  steps  requested. 

MATERIAL 

The  MATERIAL  card  is  used  to  specify  the  materials  and  physical 
parameters  to  be  used  for  the  simulation.  Material  NUMBER  1  corresponding 
to  region  number  1  is  OXIDE.  Other  insulator  specifications  allowed  are  S102, 
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NITRIDE,  SI3N4,  SAPPHIRE,  and  INSULATOr.  The  relative  permittivity 
is  appropriately  set  for  each  of  the  insulator  parameters  except  INSULATOr 
which  requires  an  explicit  PERMITTIvity  value.  It  is  optional  for  the  other 
specificatons. 

The  second  MATERIAL  card  assigns  SILICON  physical  parameters  to 
region  NUMBER  2.  Other  semiconductor  specifications  allowed  are  gallium 
arsenide  (GAAS)  and  SEMICONDuctor.  A  variety  of  physical  parameters 
are  set  by  the  SILICON  or  GAAS  parameters  or  may  be  optionally  set,  but 
they  must  be  explicitly  set  for  the  SEMICONDuctor  parameter.  N1300  and 
EG300  are  the  intrinsic  carrier  concentration  and  energy  gap  at  30(T7\,  the 
PERMITTIlivity  and  electron  AFFINITY  may  be  specified,  and  TP  and  TN 
arc  the  hole  and  electron  minority  carrier  lifetimes.  A  constant  MOBILITY 
and  a  carrier  saturation  velocity  (VSAT)  may  also  be  specified.  EGALPIIA 
and  EGBETA  are  terms  in  the  expression  for  energy  gap  variations  with 
temperature: 


Eg{T)  = 


(7Y300M 

7’/300  +  ft  ) 


(4.01) 


where  Eg  is  the  energy  gap,  T  is  the  temperature  in  °K,  a  is  EGALPIIA, 
and  0  is  EGBETA.  These  values  are  related  to  those  oT  Sze  [A.l]  by  a  = 
300a(Sze)  and  ft  =  /J(Sze)/300. 


CONTACT 

The  CONTACT  card  specifies  the  type  of  material  used  for  the  device 
electrodes.  The  NUMBER  corresponds  to  the  ELECTRODe  number. 
Alternatively,  AIL  may  be  used  to  specify  with  one  card  that  all  of  the 
contacts  use  the  same  material.  In  the  example,  the  gate  is  N.POLYSIlicon 
and  the  source  and  drain  are  ALUMINUM.  Other  allowed  materials 
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arc  P.POLYSIlicon,  MOLYBDENum,  and  TUNGSTEN,  molybdenum  dis- 
ilicide  (MO.DIS1L)  and  tungsten  disilicidc  (TU.D1SIL).  Alternatively,  the 
WORKFUNCtion  may  be  provided  explicitly.  It  is  very  useful  in  those  situa¬ 
tions  where  the  contact  characteristics  do  not  influence  device  operation  to 
specify  a  NEUTRAL  contact.  This  specification  guarantees  that  there  will  be 
no  carrier  accumulation  or  depletion  at  the  contact.  The  substrate  contact 
in  the  MOSFET  simulation,  for  example,  is  specified  in  this  manner  since  the 
simulated  substrate  contact  at  the  bottom  of  the  simulation  region  is  much 
closer  to  the  surface  than  the  actual  .substrate  contact. 

MOBILITY 

The  MOBILITY  card  is  used  to  specify  which  of  the  carrier  transport  and 
recombination  models  are  to  be  used  in  the  simulation.  This  capability  is  used 
primarily  to  aid  comparison  of  PISCES  to  other  device  simulation  programs 
which  do  not  have  the  models.  For  normal  simulations  one  would  turn  on  all 
of  the  models.  The  example  specifies  that  velocity  saturation  (YSATURAT) 
and  impurity  concentration  dependent  mobility  (CONMOB)  be  used  in  the 
simulation.  The  other  model  parameter  allowed  is  SE1JREOOM  for  Shockley- 
Read-IIall  recombination.  Each  of  these  specifications  remains  in  force  until 
terminated  with  a  NOVSATUR,  NOCONMOB,  or  NOSRHREO  specification. 

The  MOBILITY  card  is  the  last  card  in  the  SETUP  sequence.  Reading 
of  the  next  card  causes  the  initial  guess  to  be  generated  and  all  parameters 
to  be  set  as  specified. 

SOLVE. 

The  SOLVE  card  controls  the  method  of  solution  used  in  the  simulation. 
The  first.  SOLVE  card  in  the  example  specifies  that  iteration  information  (bias, 
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charge  and  current  at  each  electrode)  be  PRINTed  and  that  the  solution  be 
saved  in  a  file  called  EXOUTO  (OUTFILE).  The  PRINT  parameter  can  be 
terminated  using  NOPRINT. 

The  fourth  SOLVE  card  specifies  that  the  SJNGLEPoisson  iteration 
method  be  used  since  significant  drain  bias  is  being  applied.  The  default 
condition  is  the  MULTIPOisson  iteration  method.  One  may  also  specify 
POISSON  only  iterations  in  which  no  continuity  equation  solutions  are  per¬ 
formed.  The  default  condition  is  BOTH  where  both  sets  of  equations  are 
solved.  The  ACCELERAtion  parameter  specifies  that  the  linearizat  ion  factor 
reduction  method  of  convergence  acceleration  discussed  in  the  text  is  to  be 
used.  The  default  is  NOACCELEration.  The  alternate  acceleration  method, 
overRELAXation  is  also  available  with  the  default  ui  NORELAX.  There  are 
four  levels  of  convergence  available,  COARSE,  MEDIUM,  FINE  and  LIMIT 
with  MEDIUM  being  the  default.  Each  of  the  first  three  have  succeedingly 
tighter  convergence  limits.  The  fourth  level,  LIMIT,  specifies  that  a  number 
of  iterations  equal  to  ITLIMIT  be  executed  regardless  of  the  level  of  conver¬ 
gence.  Alternatively,  the  actual  iteration  convergence  tolerances  themselves 
P.TOLERancc  and  C.TOLERance  may  be  specified. 

The  fifth  SOLVE  card  demonstrates  its  use  in  the  bias  stepping  mode. 
No  special  consideration  is  required;  however,  the  file  name  specified  by  the 
OUTFILE  parameter  will  be  incremented  by  one  character/digit  for  each 
bias  step.  Thus  3  different  solution  files  will  be  saved  for  this  bias  stepping 
sequence;  EXOUT3,  EXOUT4  and  EXOUT5. 

The  PLOT. 2D  card  has  been  covered  earlier  with  the  card  sequence  of 
Figure  A.  1. 
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CONTOUR 

The  CONTOUR  card  is  used  to  plot  two-dimensional  contours  of  various 
device  parameters.  EquiPOTENTI A1  contours  are  plotted  at  potential 
values  of  MIN.VALUe  to  MAX.VALUe  with  DEL.VALUe  steps  Other  device 
contours  which  may  be  plotted  are  quasi-Fermi  potential  (QF.POTEN), 
DOPING,  ELECTRON  concentration,  HOLE  concentration,  net  charge 
concentration  (NET. Cl  IRC)  and  net  carrier  concentration  (NET. CARR). 
Logarithmic  contour  intervals  may  he  specified  by  LOGARITlIm.  The  values 
for  MIN.VALUe,  MAX.VALUe,  and  DEL.VALUe  are  then  the  logarithms  of 
the  desired  values.  In  order  to  plot  logarithmic  intervals  of  negative  values  of 
NET.CHRG  or  NET. CARR,  the  NEGATIVE  parameter  must  be  specified. 
The  line  type  to  be  used  in  the  contour  plot  may  be  specified  using  the 
LINE.TYPe  parameter  and  an  integer  value  between  one  and  II.  Figures  A. 4 
and  A. 5  show  the  potential  and  quasi-Fermi  potential  contour  plots  generated 
by  the  CONTOUR  cards. 
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